None Ch2_code

Table of Contents

In [ ]:
# install.packages("bizdays")
# install.packages("RQuantLib")
# install.packages("fUnitRoots")
# install.packages("TSA")
# install.packages("devtools")
# install.packages("Dict")
In [685]:
# install.packages("devtools")
also installing the dependencies ‘askpass’, ‘credentials’, ‘openssl’, ‘sys’, ‘zip’, ‘gitcreds’, ‘httr2’, ‘ini’, ‘httpuv’, ‘mime’, ‘xtable’, ‘fontawesome’, ‘sourcetools’, ‘later’, ‘promises’, ‘jquerylib’, ‘sass’, ‘systemfonts’, ‘textshaping’, ‘tinytex’, ‘diffobj’, ‘rematch2’, ‘clipr’, ‘gert’, ‘gh’, ‘purrr’, ‘rappdirs’, ‘rprojroot’, ‘rstudioapi’, ‘whisker’, ‘cachem’, ‘shiny’, ‘callr’, ‘processx’, ‘bslib’, ‘downlit’, ‘httr’, ‘ragg’, ‘rmarkdown’, ‘xml2’, ‘htmlwidgets’, ‘stringr’, ‘prettyunits’, ‘xopen’, ‘brew’, ‘commonmark’, ‘stringi’, ‘cpp11’, ‘brio’, ‘praise’, ‘ps’, ‘waldo’, ‘usethis’, ‘desc’, ‘fs’, ‘memoise’, ‘miniUI’, ‘pkgbuild’, ‘pkgdown’, ‘pkgload’, ‘profvis’, ‘rcmdcheck’, ‘remotes’, ‘roxygen2’, ‘rversions’, ‘sessioninfo’, ‘testthat’, ‘urlchecker’


Warning message in install.packages("devtools"):
“installation of package ‘gert’ had non-zero exit status”
Warning message in install.packages("devtools"):
“installation of package ‘usethis’ had non-zero exit status”
Warning message in install.packages("devtools"):
“installation of package ‘devtools’ had non-zero exit status”
Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

In [182]:
# install.packages("Dict")
also installing the dependencies ‘tidyselect’, ‘dplyr’


Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

In [15]:
# devtools::create("AFTSCode")
In [1]:
library(fUnitRoots)
library(bizdays, RQuantLib)
library(TSA) # eacf
library(forecast)
print(packageVersion("bizdays"))
print(packageVersion("RQuantLib"))
library(knitr)
library(IRdisplay)
library(fBasics)
library(Dict)
library(sandwich) # Use heteroskedasticity-consistent (HC) variance estimation
library(lmtest) # Apply HC variance to linear models
library(tseries) # Jarque-Bera test
Attaching package: ‘bizdays’


The following object is masked from ‘package:stats’:

    offset



Attaching package: ‘TSA’


The following objects are masked from ‘package:stats’:

    acf, arima


The following object is masked from ‘package:utils’:

    tar


Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Registered S3 methods overwritten by 'forecast':
  method       from
  fitted.Arima TSA 
  plot.Arima   TSA 


Attaching package: ‘forecast’


The following object is masked from ‘package:bizdays’:

    bizdays


[1] ‘1.0.16’
[1] ‘0.4.21’
Attaching package: ‘fBasics’


The following objects are masked from ‘package:TSA’:

    kurtosis, skewness


Loading required package: zoo


Attaching package: ‘zoo’


The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric


In [13]:
# detach("package:AFTSCode", unload = TRUE) # First detach the package
# unloadNamespace("AFTSCode")
In [14]:
library(AFTSCode)
In [19]:
# find("plot_pacf_acf"); find("draw_arima_forecast_fig")
  1. '.GlobalEnv'
  2. 'package:AFTSCode'
'package:AFTSCode'
In [1]:
# install_redir_output <- function(lib, log_fname) {
#     # Open a file connection for writing
#     sink(log_fname, append = TRUE)

#     # Install packages (output will be redirected to the file)
#     install.packages(lib)

#     # Close the file connection
#     sink()
# }

# plot_time_fig <- function(ts, main, xlab, ylab) {
#     par(bg = "white")
#     plot(ts, type = 'l', main = main, xlab = xlab, ylab = ylab)
#     points(ts, pch = '*')
# }

# cal_phi_0 <- function(arima_md, ord, digits=6) {
#     phi_0 = as.numeric((1-sum(arima_md$coef[1:ord]))*arima_md$coef['intercept'])
#     print(format(phi_0, digits = digits))
#     phi_0
# }

# get_mu <- function(arima_mod, digits=6) {
#     mu = as.numeric(arima_mod$coef['intercept'])
#     print(format(mu, digits = digits))
#     mu
# }

# plot_arima_forecast_fig <- function(
#     da_ts, eotr, h, npts, frequency, 
#     order, fixed, method, transform.pars,
#     main, xlab, ylab, ylim = NULL
# ) {
#     par(bg = "white")
#     # arima model
#     tr_da_ts = ts(da_ts[1:eotr], frequency = frequency, start = start(da_ts))
#     if (is.null(transform.pars)) {
#         ts_fm3 = arima(tr_da_ts, order = order, fixed = fixed, method = method)
#     } else {
#         ts_fm3 = arima(tr_da_ts, order = order, fixed = fixed, method = method, transform.pars = transform.pars)
#     }
#     print(ts_fm3$nobs)
#     # Forecast
#     ts_fm3$x <- tr_da_ts # https://stackoverflow.com/a/42464130/4307919
#     ts_fc_res = forecast(ts_fm3, h = h)
#     # Plot forecast
#     xmin = time(da_ts)[eotr]-npts/frequency
# #     xmax = end(da_ts)[1]+(end(da_ts)[2]+2)/frequency
#     xmax = time(da_ts)[eotr]+(max(h, length(da_ts)-eotr)+1)/frequency
#     cat(xmin, ";", xmax)
    
#     plot(ts_fc_res, xlim = c(xmin, xmax), ylim = ylim, main = main, xlab = xlab, ylab = ylab)
#     # Plot forecast mean
#     dummy_1st_fmean_ts = ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$mean)), frequency = frequency, start = end(tr_da_ts))
#     lines(dummy_1st_fmean_ts)
#     points(dummy_1st_fmean_ts, pch = 1)
#     # Plot confidence interval
#     dummy_1st_flower_ts = ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$lower[,2])), frequency = frequency, start = end(tr_da_ts))
#     dummy_1st_fupper_ts = ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$upper[,2])), frequency = frequency, start = end(tr_da_ts))
#     lines(dummy_1st_flower_ts, lty=2)
#     lines(dummy_1st_fupper_ts, lty=2)
#     # Plot original data
#     orig_plot_ts = ts(da_ts[(eotr-npts+1):length(da_ts)], frequency = frequency, start = time(da_ts)[eotr]-(npts-1)/frequency)
#     lines(orig_plot_ts)
#     points(orig_plot_ts, pch = 19)
#     ts_fc_res
# }

# comb_forecast_res <- function(forecast_obj, da_ts, eotr, freq) {
#     display(summary(forecast_obj))
#     fc_std_err = (forecast_obj$upper[,2]-forecast_obj$lower[,2])/2/qnorm(p = 0.975)
#     actual_ts = ts(da_ts[(eotr+1):length(da_ts)], frequency = freq, start = time(da_ts)[eotr+1])
#     display(forecast_obj$mean); display(fc_std_err); display(actual_ts)
#     multistep_ahead_forecast_tb = cbind(forecast_obj$mean, fc_std_err, actual_ts)
#     dimnames(multistep_ahead_forecast_tb)[[2]] <- c("Forecast", "Std. Error", "Actual")
#     multistep_ahead_forecast_tb
# }

# plot_unit_root_figs <- function(da_ts, freq, start, lag_max, main, xlab, ylab) {
#     par(bg = "white")
# #     lgdp_ts = ts(lgdp, frequency = 4, start = c(1947 ,1))
#     diff_da_ts = ts(diff(da_ts), frequency = freq, start = start)
#     plot(da_ts, col = 'black', type = 'l', main = main, xlab = xlab, ylab = ylab)
#     acf(da_ts, lag.max = 20)
#     plot_time_fig(diff_da_ts, main = main, xlab = xlab, ylab = paste("diff(", ylab, ")", sep = ""))
#     pacf(diff_da_ts, lag.max = 20)
# }

# plot_acf <- function(da, lag.max = NULL, main=NULL, w=NULL, h=NULL, ...) {
#     if (is.null(w) | is.null(h)) {
#         par(bg = 'white')
#     } else {
#         par(bg = 'white', pin = c(w, h))
#     }
#     plot(acf(da, lag.max = lag.max, plot = F, `drop.lag.0` = F, ...), main = main)
# }

# plot_pacf_acf <- function(da, lag.max = NULL, main=NULL, w=NULL, h=NULL, ...) {
#     if (is.null(w) | is.null(h)) {
#         par(mfrow = c(2, 1), bg = 'white')
#     } else {
#         par(mfrow = c(2, 1), bg = 'white', pin = c(w, h))
#     }
#     plot(pacf(da, lag.max = lag.max, plot = F, `drop.lag.0` = F, ...), main = main)
#     plot(acf(da, lag.max = lag.max, plot = F, `drop.lag.0` = F, ...), main = main)
# }
In [ ]:

Sec. 2.2

Box-Ljung test (IBM and CRSP Value-Weighted Index)

  • IBM monthly returns shows almost no serial correlation.
  • CRSP Value-Weighted Index monthly returns shows medium momentum at lag-1 and medium reverting at lag-3.
In [3]:
da = read.table("../AFTS_data/Ch02/m-ibm3dx2608.txt", header=T)
da[1:5,]
A data.frame: 5 × 5
dateibmrtnvwrtnewrtnsprtn
<int><dbl><dbl><dbl><dbl>
119260130-0.010381 0.000724 0.023174 0.022472
219260227-0.024476-0.033374-0.053510-0.043956
319260331-0.115591-0.064341-0.096824-0.059113
419260430 0.089783 0.038358 0.032946 0.022688
519260528 0.036932 0.012172 0.001035 0.007679
In [9]:
find("acf")
  1. 'package:TSA'
  2. 'package:stats'
In [22]:
# IBM (Fig. 2.1)
plot_pacf_acf(da$ibmrtn, lag.max = 100)
No description has been provided for this image
No description has been provided for this image
In [24]:
# CRSP Value-Weighted Index (Fig. 2.2)
plot_pacf_acf(da$vwrtn, lag.max = 100)
No description has been provided for this image
No description has been provided for this image
In [59]:
Box.test(da$ibmrtn, lag = 5, type = 'Ljung')
	Box-Ljung test

data:  da$ibmrtn
X-squared = 3.3682, df = 5, p-value = 0.6434
In [60]:
libm = log(da$ibmrtn+1)
Box.test(libm, lag = 5, type = 'Ljung')
	Box-Ljung test

data:  libm
X-squared = 3.5236, df = 5, p-value = 0.6198

Find IBM monthly and yearly expected return

  • For the log(Return)
    • $\mu_{monthly}\approx 0.010891$
    • $\sigma_{monthly}\approx 0.070333$
In [81]:
# µ, sigma, and monthly price ratio
ibm_basic_stats <- apply(as.data.frame(libm), 2, basicStats)
ibm_basic_stats$libm['Mean',];
ibm_basic_stats$libm['Stdev',];
mnthly_ibm_pr_ratio = exp(ibm_basic_stats$libm['Mean',]+ibm_basic_stats$libm['Stdev',]^2/2)
mnthly_ibm_pr_ratio; mnthly_ibm_pr_ratio^12;
0.010891
0.070333
1.01345406773567
1.17394795060543

Sec. 2.4

AR(3)

In [9]:
gnp = scan(file = "../AFTS_data/Ch02/dgnp82.txt")
gnp1 = ts(gnp, frequency = 4, start = c(1947, 2))
In [25]:
plot_time_fig(gnp1, main = "Gross National Product", xlab = "Year", ylab = "Growth")
No description has been provided for this image
  • Find the AR order using AIC
In [10]:
m1 = ar(gnp, method = 'mle')
m1$order
3
In [11]:
m2 = arima(gnp, order = c(3, 0, 0))
m2
Call:
arima(x = gnp, order = c(3, 0, 0))

Coefficients:
         ar1     ar2      ar3  intercept
      0.3480  0.1793  -0.1423     0.0077
s.e.  0.0745  0.0778   0.0745     0.0012

sigma^2 estimated as 9.427e-05:  log likelihood = 565.84,  aic = -1121.68
  • $\phi_0 = \mu*(1-\phi_1-\phi_2-\phi_3)$
In [15]:
cal_phi_0(m2, ord=3)
[1] "0.00472311"
0.00472311193618254

Residual standar error

In [52]:
sqrt(m2$sigma2)
0.00970932215017931
  • Characteristic equation
In [57]:
?polyroot
In [53]:
p1 = c(1, -m2$coef[1:3])
roots = polyroot(p1)
roots
  1. 1.59025281352281+1.06388168563059i
  2. -1.92015204104451-0i
  3. 1.59025281352281-1.06388168563059i
In [71]:
for (i in 1:3) {
    print(i)
    print(as.numeric(Mod(1-m2$coef[1]*roots[i]-m2$coef[2]*roots[i]^2-m2$coef[3]*roots[i]^3)))
}
[1] 1
[1] 2.237726e-16
[1] 2
[1] 4.641385e-16
[1] 3
[1] 2.237726e-16
  • Modulus of the roots
In [77]:
roots; Mod(roots)
  1. 1.59025281352281+1.06388168563059i
  2. -1.92015204104451-0i
  3. 1.59025281352281-1.06388168563059i
  1. 1.91330819575347
  2. 1.92015204104451
  3. 1.91330819575347
In [76]:
1/roots; Mod(1/roots)
  1. 0.434406493995099-0.290618641985979i
  2. -0.520792092826164+0i
  3. 0.434406493995099+0.290618641985979i
  1. 0.522654950320849
  2. 0.520792092826164
  3. 0.522654950320849
  • Compute the average length of business cycles
  • Doesn't matter if using the characteristic roots or the inverse of them.
  • Unit is "quarter". So average business cycle is about 10.66 quarters.
In [80]:
k = 2*pi/acos(0.434406493995099/0.522654950320849)
k
10.6563777674414
In [81]:
k = 2*pi/acos(1.59025281352281/1.91330819575347)
k
10.6563777674414

AIC

In [84]:
gnp = scan(file = '../AFTS_data/Ch02/dgnp82.txt')
ord = ar(gnp, method = 'mle')
format(ord$aic, digits = 6)
0
'27.846690'
1
' 2.741632'
2
' 1.603242'
3
' 0.000000'
4
' 0.302785'
5
' 2.242661'
6
' 4.052084'
7
' 6.025475'
8
' 5.904668'
9
' 7.571863'
10
' 7.895334'
11
' 9.678873'
12
' 7.197545'
In [86]:
ord$order; ord$order.max
3
12

AR(3) example

In [66]:
da = read.table("../AFTS_data/Ch02/m-ibm3dx2608.txt", header = T)
da[1:5,]
A data.frame: 5 × 5
dateibmrtnvwrtnewrtnsprtn
<int><dbl><dbl><dbl><dbl>
119260130-0.010381 0.000724 0.023174 0.022472
219260227-0.024476-0.033374-0.053510-0.043956
319260331-0.115591-0.064341-0.096824-0.059113
419260430 0.089783 0.038358 0.032946 0.022688
519260528 0.036932 0.012172 0.001035 0.007679

Parameter estimation for value-weighted index

  • Try simple return
In [26]:
vwrtn_m3 = arima(da$vwrtn, order = c(3, 0, 0))
vwrtn_m3
Call:
arima(x = da$vwrtn, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1158  -0.0187  -0.1042     0.0089
s.e.  0.0315   0.0317   0.0317     0.0017

sigma^2 estimated as 0.002875:  log likelihood = 1500.86,  aic = -2991.73
In [67]:
vwrtn_ts = ts(da$vwrtn, frequency = 12, start = c(1926, 1))
plot_time_fig(vwrtn_ts, main = "Value-weighted index return", xlab = "Year", ylab = "Return")
No description has been provided for this image
In [63]:
est_mrtn = get_mu(vwrtn_m3)
(1+est_mrtn)^12-1
[1] "0.00894788"
0.112819646650081
  • Try log-return
In [69]:
lvwrtn_ts = ts(log(1+da$vwrtn), frequency = 12, start = c(1926, 1))
lvwrtn_m3 = arima(lvwrtn_ts, order = c(3, 0, 0))
lvwrtn_m3
Call:
arima(x = lvwrtn_ts, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1105  -0.0139  -0.0953     0.0074
s.e.  0.0315   0.0317   0.0318     0.0017

sigma^2 estimated as 0.002886:  log likelihood = 1498.89,  aic = -2989.78
In [70]:
plot_time_fig(lvwrtn_ts, main = "Value-weighted index log-return", xlab = "Year", ylab = "log(Return)")
No description has been provided for this image
In [65]:
est_mlrtn = get_mu(lvwrtn_m3)
exp(est_mlrtn*12)-1
[1] "0.00744436"
0.0934439038357271

Calculate the actual yearly return

  • The 9.3% per annum is close to 11.3% per annum from the estimated model on simple return.
  • If I transfer the simple return to log return using log(1+rtn), the estimated return per annum is 9.3%, which is closer to the actual value.
In [90]:
vw = da$vwrtn # Value-weighted index
t1 = prod(vw+1) # Simple gross return
t1; length(vw); t1^(12/length(vw))-1; # monthly return per annum
1592.95348104506
996
0.0929008410097438

Model checking

In [40]:
m3 = arima(da$vwrtn, order = c(3, 0, 0))
m3
Call:
arima(x = da$vwrtn, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1158  -0.0187  -0.1042     0.0089
s.e.  0.0315   0.0317   0.0317     0.0017

sigma^2 estimated as 0.002875:  log likelihood = 1500.86,  aic = -2991.73
  • Compute the intercept $\phi_0 = \mu*(1-\phi_1-\phi_2-\phi_3)$
In [41]:
cal_phi_0(m3, ord=3)
[1] "0.00901194"
0.00901193749698098
  • Compute the standard error of residuals
In [42]:
attributes(m3)
$names
  1. 'coef'
  2. 'sigma2'
  3. 'var.coef'
  4. 'mask'
  5. 'loglik'
  6. 'aic'
  7. 'arma'
  8. 'residuals'
  9. 'call'
  10. 'series'
  11. 'code'
  12. 'n.cond'
  13. 'nobs'
  14. 'model'
$class
'Arima'
In [43]:
c(
    sqrt(sum(m3$residuals^2)/(m3$nobs-2*length(m3$coef)-1)),
    sqrt(sum(m3$residuals^2)/(m3$nobs-length(m3$coef))),
    sqrt(sum(m3$residuals^2)/m3$nobs)
)
  1. 0.0538628065685737
  2. 0.0537268921289658
  3. 0.0536188982667509
In [44]:
sqrt(m3$sigma2)
0.0536188982667509
  • Ljung-Box test (R uses 12 dof)
In [45]:
Box.test(m3$residuals, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  m3$residuals
X-squared = 16.352, df = 12, p-value = 0.1756
  • Compute p-value using 9 dof (12-3=9, the parameter for AR does not include the intercept.)
In [122]:
pv = 1-pchisq(16.35,9)
pv
0.0599227642192111
  • Refine the model to rule out statistically insignificant coef by fixing the AR(2) coef to zero (pp51)
In [362]:
mm3 = arima(da$vwrtn, order = c(3,0,0), fixed = c(NA,0,NA,NA))
mm3
Warning message in arima(da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA)):
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA))

Coefficients:
         ar1  ar2      ar3  intercept
      0.1136    0  -0.1063     0.0089
s.e.  0.0313    0   0.0315     0.0017

sigma^2 estimated as 0.002876:  log likelihood = 1500.69,  aic = -2993.38
In [366]:
mm3_1 = arima(da$vwrtn, order = c(3,0,0), fixed = c(NA,0,NA,NA), method = "ML")
mm3_1; cal_phi_0(mm3_1, ord = 3)
Warning message in arima(da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), method = "ML")

Coefficients:
         ar1  ar2      ar3  intercept
      0.1136    0  -0.1063     0.0089
s.e.  0.0313    0   0.0315     0.0017

sigma^2 estimated as 0.002876:  log likelihood = 1500.69,  aic = -2993.38
[1] "0.0088825"
0.00888250044781528
In [367]:
mm3_2 = arima(da$vwrtn, order = c(3,0,0), fixed = c(NA,0,NA,NA), method = "CSS")
mm3_2; cal_phi_0(mm3_2, ord = 3)
Warning message in arima(da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), method = "CSS")

Coefficients:
         ar1  ar2      ar3  intercept
      0.1126    0  -0.1064     0.0091
s.e.  0.0313    0   0.0315     0.0017

sigma^2 estimated as 0.002878:  part log likelihood = 1500.33
[1] "0.00899938"
0.00899938186119909
In [368]:
mm3_3 = arima(da$vwrtn, order = c(3,0,0), fixed = c(NA,0,NA,NA), method = "CSS-ML")
mm3_3; cal_phi_0(mm3_3, ord = 3)
Warning message in arima(da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), method = "CSS-ML")

Coefficients:
         ar1  ar2      ar3  intercept
      0.1136    0  -0.1063     0.0089
s.e.  0.0313    0   0.0315     0.0017

sigma^2 estimated as 0.002876:  log likelihood = 1500.69,  aic = -2993.38
[1] "0.00888153"
0.00888153260242181
In [47]:
cal_phi_0(mm3, ord = 3)
[1] "0.00888153"
0.00888153260242181
  • In the arima software, when estimating the residual variance, they use sample size $T$, instead of $T-2*p-1$.
In [56]:
sum(mm3$residuals^2)/(length(da$vwrtn)); mm3$sigma2; sqrt(mm3$sigma2)
0.00287599654600986
0.00287599654600986
0.0536283185081339
In [49]:
box_test <- Box.test(mm3$residuals, lag = 12, type = 'Ljung')
box_test
	Box-Ljung test

data:  mm3$residuals
X-squared = 16.828, df = 12, p-value = 0.1562
In [50]:
attributes(box_test)
$names
  1. 'statistic'
  2. 'parameter'
  3. 'p.value'
  4. 'method'
  5. 'data.name'
$class
'htest'
  • Compute p-value using 10 dof (12-2=10, the parameter for AR does not include the intercept.)
In [51]:
pv = as.numeric(1-pchisq(box_test$statistic, 10))
pv; box_test$statistic
0.0782661025605844
X-squared: 16.8276328679107

AR(3) forecasting example

  • See issues related to scoping for forecast: link
In [70]:
# install.packages("forecast")
also installing the dependencies ‘farver’, ‘labeling’, ‘munsell’, ‘R6’, ‘RColorBrewer’, ‘viridisLite’, ‘pkgconfig’, ‘xts’, ‘TTR’, ‘curl’, ‘gtable’, ‘isoband’, ‘scales’, ‘tibble’, ‘withr’, ‘quadprog’, ‘quantmod’, ‘colorspace’, ‘fracdiff’, ‘generics’, ‘ggplot2’, ‘lmtest’, ‘magrittr’, ‘tseries’, ‘RcppArmadillo’


Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

In [644]:
da = read.table("../AFTS_data/Ch02/m-ibm3dx2608.txt", header = T)
da[1:5,]
A data.frame: 5 × 5
dateibmrtnvwrtnewrtnsprtn
<int><dbl><dbl><dbl><dbl>
119260130-0.010381 0.000724 0.023174 0.022472
219260227-0.024476-0.033374-0.053510-0.043956
319260331-0.115591-0.064341-0.096824-0.059113
419260430 0.089783 0.038358 0.032946 0.022688
519260528 0.036932 0.012172 0.001035 0.007679
In [651]:
f_da_ts = ts(da$vwrtn[1:tr_n], frequency = 12, start = c(1926, 1))
ts_fm3 = arima(f_da_ts, order = c(3,0,0))
# ts_fc_res = forecast(ts_fm3, h = 12)
ts_fc_res = predict(ts_fm3, n.ahead=12)
ts_fc_res
$pred
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20080.0074556860.0159524880.0117016630.0098066530.0087708860.0091646440.0094325200.0095650950.0095305550.0094951510.0094777480.009480420
$se
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20080.053289020.053573400.053575740.053906750.053920260.053920290.053924560.053924930.053924930.053924990.053925000.05392500
In [660]:
f_da_ts = ts(da$vwrtn[1:tr_n], frequency = 12, start = c(1926, 1))
ts_fm3 = arima(f_da_ts, order = c(3,0,0))
ts_fm3$x <- f_da_ts
ts_fc_res = forecast(ts_fm3, h = 12)
summary(ts_fc_res)
Forecast method: ARIMA(3,0,0) with non-zero mean

Model Information:

Call:
arima(x = f_da_ts, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1034  -0.0201  -0.1089     0.0095
s.e.  0.0317   0.0318   0.0317     0.0017

sigma^2 estimated as 0.00284:  log likelihood = 1488.85,  aic = -2969.71

Error measures:
                        ME       RMSE        MAE      MPE     MAPE      MASE
Training set -1.234691e-05 0.05328902 0.03757996 101.6395 166.5597 0.7001715
                    ACF1
Training set 0.002856914

Forecasts:
         Point Forecast       Lo 80      Hi 80       Lo 95     Hi 95
Jan 2008    0.007455686 -0.06083694 0.07574831 -0.09698887 0.1119002
Feb 2008    0.015952488 -0.05270459 0.08460957 -0.08904946 0.1209544
Mar 2008    0.011701663 -0.05695841 0.08036174 -0.09330486 0.1167082
Apr 2008    0.009806653 -0.05927763 0.07889093 -0.09584863 0.1154619
May 2008    0.008770886 -0.06033071 0.07787249 -0.09691089 0.1144527
Jun 2008    0.009164644 -0.05993698 0.07826627 -0.09651718 0.1148465
Jul 2008    0.009432520 -0.05967458 0.07853962 -0.09625768 0.1151227
Aug 2008    0.009565095 -0.05954249 0.07867268 -0.09612583 0.1152560
Sep 2008    0.009530555 -0.05957703 0.07863814 -0.09616037 0.1152215
Oct 2008    0.009495151 -0.05961250 0.07860281 -0.09619589 0.1151862
Nov 2008    0.009477748 -0.05962992 0.07858541 -0.09621331 0.1151688
Dec 2008    0.009480420 -0.05962725 0.07858808 -0.09621063 0.1151715
In [652]:
cal_phi_0(ts_fm3, ord = 3)
[1] "0.00972842"
0.00972842337769138
In [653]:
sqrt(ts_fm3$sigma2)
0.0532890177490275
In [661]:
# Use forecast
par(bg = "white")
xend = end(f_da_ts)[1]+end(f_da_ts)[2]/12
xmin = xend-1
xmax = xend+1
plot(ts_fc_res, xlim = c(xmin, xmax))
No description has been provided for this image
In [657]:
# Use predict
par(bg = "white")
xend = end(f_da_ts)[1]+end(f_da_ts)[2]/12
xmin = xend-1
xmax = xend+1
plot(ts_fc_res$pred, xlim = c(xmin, xmax), ylim = c(-0.3, 0.4))
No description has been provided for this image
In [281]:
attributes(ts_fc_res)
$names
  1. 'method'
  2. 'model'
  3. 'level'
  4. 'mean'
  5. 'lower'
  6. 'upper'
  7. 'x'
  8. 'series'
  9. 'fitted'
  10. 'residuals'
$class
'forecast'
  • Make multi-steps forecast plots (Fig. 2.7)
In [421]:
freq = 12
vwrtn_ts = ts(da$vwrtn, frequency = freq, start = c(1926, 1))
eotr=984
ts_fc_res = plot_arima_forecast_fig(
    da_ts=vwrtn_ts, 
    eotr=eotr, h=12, npts=12, frequency=freq, 
    order=c(3,0,0), fixed=NULL, method=NULL, transform.pars=TRUE,
    main="Forecasts from ARIMA(3,0,0) with non-zero mean for ewrtn (Fig. 2.7)", 
    xlab="Year", ylab="Return", ylim=c(-0.21, 0.21)
)
2006.917 ; 2009
No description has been provided for this image
  • Forecast values (Table 2.2)
In [292]:
summary(ts_fc_res)
Forecast method: ARIMA(3,0,0) with non-zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1034  -0.0201  -0.1089     0.0095
s.e.  0.0317   0.0318   0.0317     0.0017

sigma^2 estimated as 0.00284:  log likelihood = 1488.85,  aic = -2967.71

Error measures:
                        ME       RMSE        MAE      MPE     MAPE      MASE
Training set -1.234691e-05 0.05328902 0.03757996 101.6395 166.5597 0.7001715
                    ACF1
Training set 0.002856914

Forecasts:
         Point Forecast       Lo 80      Hi 80       Lo 95     Hi 95
Jan 2008    0.007455686 -0.06083694 0.07574831 -0.09698887 0.1119002
Feb 2008    0.015952488 -0.05270459 0.08460957 -0.08904946 0.1209544
Mar 2008    0.011701663 -0.05695841 0.08036174 -0.09330486 0.1167082
Apr 2008    0.009806653 -0.05927763 0.07889093 -0.09584863 0.1154619
May 2008    0.008770886 -0.06033071 0.07787249 -0.09691089 0.1144527
Jun 2008    0.009164644 -0.05993698 0.07826627 -0.09651718 0.1148465
Jul 2008    0.009432520 -0.05967458 0.07853962 -0.09625768 0.1151227
Aug 2008    0.009565095 -0.05954249 0.07867268 -0.09612583 0.1152560
Sep 2008    0.009530555 -0.05957703 0.07863814 -0.09616037 0.1152215
Oct 2008    0.009495151 -0.05961250 0.07860281 -0.09619589 0.1151862
Nov 2008    0.009477748 -0.05962992 0.07858541 -0.09621331 0.1151688
Dec 2008    0.009480420 -0.05962725 0.07858808 -0.09621063 0.1151715
In [445]:
install_redir_output(lib = "knitr", log_fname = "knitr_install.log")
also installing the dependencies ‘highr’, ‘xfun’, ‘yaml’


Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

In [502]:
vwrtn_multistep_ahead_forecast_tb <- comb_forecast_res(ts_fc_res, vwrtn_ts, eotr, freq)
vwrtn_multistep_ahead_forecast_tb
Forecast method: ARIMA(3,0,0) with non-zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, fixed = fixed, method = method)

Coefficients:
         ar1      ar2      ar3  intercept
      0.1034  -0.0201  -0.1089     0.0095
s.e.  0.0317   0.0318   0.0317     0.0017

sigma^2 estimated as 0.00284:  log likelihood = 1488.85,  aic = -2967.71

Error measures:
                        ME       RMSE        MAE      MPE     MAPE      MASE
Training set -1.234691e-05 0.05328902 0.03757996 101.6395 166.5597 0.7001715
                    ACF1
Training set 0.002856914

Forecasts:
         Point Forecast       Lo 80      Hi 80       Lo 95     Hi 95
Jan 2008    0.007455686 -0.06083694 0.07574831 -0.09698887 0.1119002
Feb 2008    0.015952488 -0.05270459 0.08460957 -0.08904946 0.1209544
Mar 2008    0.011701663 -0.05695841 0.08036174 -0.09330486 0.1167082
Apr 2008    0.009806653 -0.05927763 0.07889093 -0.09584863 0.1154619
May 2008    0.008770886 -0.06033071 0.07787249 -0.09691089 0.1144527
Jun 2008    0.009164644 -0.05993698 0.07826627 -0.09651718 0.1148465
Jul 2008    0.009432520 -0.05967458 0.07853962 -0.09625768 0.1151227
Aug 2008    0.009565095 -0.05954249 0.07867268 -0.09612583 0.1152560
Sep 2008    0.009530555 -0.05957703 0.07863814 -0.09616037 0.1152215
Oct 2008    0.009495151 -0.05961250 0.07860281 -0.09619589 0.1151862
Nov 2008    0.009477748 -0.05962992 0.07858541 -0.09621331 0.1151688
Dec 2008    0.009480420 -0.05962725 0.07858808 -0.09621063 0.1151715
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20080.0074556860.0159524880.0117016630.0098066530.0087708860.0091646440.0094325200.0095650950.0095305550.0094951510.0094777480.009480420
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20080.053289020.053573400.053575740.053906750.053920260.053920290.053924560.053924930.053924930.053924990.053925000.05392500
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
2008-0.062305-0.022040-0.010470 0.051144 0.023825-0.078625-0.013150 0.011042-0.098060-0.184726-0.085206 0.021482
A Time Series: 12 × 3
ForecastStd. ErrorActual
Jan 20080.0074556860.05328902-0.062305
Feb 20080.0159524880.05357340-0.022040
Mar 20080.0117016630.05357574-0.010470
Apr 20080.0098066530.05390675 0.051144
May 20080.0087708860.05392026 0.023825
Jun 20080.0091646440.05392029-0.078625
Jul 20080.0094325200.05392456-0.013150
Aug 20080.0095650950.05392493 0.011042
Sep 20080.0095305550.05392493-0.098060
Oct 20080.0094951510.05392499-0.184726
Nov 20080.0094777480.05392500-0.085206
Dec 20080.0094804200.05392500 0.021482
In [497]:
class(vwrtn_multistep_ahead_forecast_tb); attributes(vwrtn_multistep_ahead_forecast_tb)
  1. 'mts'
  2. 'ts'
  3. 'matrix'
  4. 'array'
$dim
  1. 12
  2. 3
$dimnames
  1. NULL
    1. 'Forecast'
    2. 'Std. Error'
    3. 'Actual'
$tsp
  1. 2008
  2. 2008.91666666667
  3. 12
$class
  1. 'mts'
  2. 'ts'
  3. 'matrix'
  4. 'array'

Find Volume-weighted monthly and yearly expected return

  • For the log(Return)
    • $\mu_{monthly}\approx 0.007403$
    • $\sigma_{monthly}\approx 0.05434$
In [74]:
# µ, sigma, monthly pr ratio, yearly pr ratio
vwrtn_basic_stats = apply(as.data.frame(lvwrtn_ts), 2, basicStats)
vwrtn_basic_stats$x['Mean',];
vwrtn_basic_stats$x['Stdev',];
mnthly_vw_pr_ratio = exp(vwrtn_basic_stats$x['Mean',]+vwrtn_basic_stats$x['Stdev',]^2/2);
mnthly_vw_pr_ratio; mnthly_vw_pr_ratio^12;
0.007403
0.05434
1.00891895677127
1.11243689916729

Sec. 2.5

  • See arima doc: link

MA(9) fitting with fixed coefs example (Sec. 2.5.3)

In [77]:
da = read.table("../AFTS_data/Ch02/m-ibm3dx2608.txt", header = T)
da[1:5,]
A data.frame: 5 × 5
dateibmrtnvwrtnewrtnsprtn
<int><dbl><dbl><dbl><dbl>
119260130-0.010381 0.000724 0.023174 0.022472
219260227-0.024476-0.033374-0.053510-0.043956
319260331-0.115591-0.064341-0.096824-0.059113
419260430 0.089783 0.038358 0.032946 0.022688
519260528 0.036932 0.012172 0.001035 0.007679
In [78]:
ewrtn_ts = ts(da$ewrtn, frequency = 12, start = c(1926, 1))
plot_time_fig(ts = ewrtn_ts, main = "Equal-weighted index return (Fig. 2.8)", xlab = "Year", ylab = "Simple Return")
No description has been provided for this image
In [416]:
par(bg = "white")
acf(da$ewrtn, main = "ACF: Equal-weighted index return (Fig. 2.8)")
No description has been provided for this image
In [360]:
find("arima")
'package:stats'
  • Fit MA(9) process
In [388]:
ewrtn_ma_mod = arima(
    da$ewrtn, 
    order = c(0,0,9), 
    fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA), # Include the mean
    transform.pars = F # Doesn't make difference
)
ewrtn_ma_mod; sqrt(ewrtn_ma_mod$sigma2)
Call:
arima(x = da$ewrtn, order = c(0, 0, 9), transform.pars = F, fixed = c(NA, 0, 
    NA, 0, 0, 0, 0, 0, NA, NA))

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1909    0  -0.1199    0    0    0    0    0  0.1227     0.0122
s.e.  0.0293    0   0.0338    0    0    0    0    0  0.0312     0.0027

sigma^2 estimated as 0.005097:  log likelihood = 1215.61,  aic = -2421.22
0.0713936744364429
In [389]:
ewrtn_ma_mod_1 = arima(
    da$ewrtn, 
    order = c(0,0,9), 
    fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA),
    method = 'ML',
    transform.pars = F # Doesn't make difference
)
ewrtn_ma_mod_1; sqrt(ewrtn_ma_mod_1$sigma2)
Call:
arima(x = da$ewrtn, order = c(0, 0, 9), transform.pars = F, fixed = c(NA, 0, 
    NA, 0, 0, 0, 0, 0, NA, NA), method = "ML")

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1909    0  -0.1199    0    0    0    0    0  0.1227     0.0122
s.e.  0.0293    0   0.0338    0    0    0    0    0  0.0312     0.0027

sigma^2 estimated as 0.005097:  log likelihood = 1215.61,  aic = -2421.22
0.0713936752812398
In [390]:
ewrtn_ma_mod_2 = arima(
    da$ewrtn, 
    order = c(0,0,9), 
    fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA),
    method = 'CSS',
    transform.pars = F # Doesn't make difference
)
ewrtn_ma_mod_2; sqrt(ewrtn_ma_mod_2$sigma2)
Call:
arima(x = da$ewrtn, order = c(0, 0, 9), transform.pars = F, fixed = c(NA, 0, 
    NA, 0, 0, 0, 0, 0, NA, NA), method = "CSS")

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1909    0  -0.1198    0    0    0    0    0  0.1232     0.0122
s.e.  0.0293    0   0.0338    0    0    0    0    0  0.0313     0.0027

sigma^2 estimated as 0.005098:  part log likelihood = 1215.65
0.0713987962817682
In [394]:
ewrtn_ma_mod_3 = arima(
    da$ewrtn, 
    order = c(0,0,9), 
    fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA),
    method = 'CSS-ML',
    transform.pars = F # Doesn't make difference
)
ewrtn_ma_mod_3; sqrt(ewrtn_ma_mod_3$sigma2)
Call:
arima(x = da$ewrtn, order = c(0, 0, 9), transform.pars = F, fixed = c(NA, 0, 
    NA, 0, 0, 0, 0, 0, NA, NA), method = "CSS-ML")

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1909    0  -0.1199    0    0    0    0    0  0.1227     0.0122
s.e.  0.0293    0   0.0338    0    0    0    0    0  0.0312     0.0027

sigma^2 estimated as 0.005097:  log likelihood = 1215.61,  aic = -2421.22
0.0713936744364429
  • Box-Ljung test (us 12 dof)
In [391]:
box_test_res = Box.test(ewrtn_ma_mod$residuals, lag = 12, type = 'Ljung');
box_test_res; 1-pchisq(box_test_res$statistic, df = 9); 1-pchisq(box_test_res$statistic, df = 12)
	Box-Ljung test

data:  ewrtn_ma_mod$residuals
X-squared = 17.604, df = 12, p-value = 0.1283
X-squared: 0.0400560859710883
X-squared: 0.128254552132693
In [392]:
box_test_res_1 = Box.test(ewrtn_ma_mod_1$residuals, lag = 12, type = 'Ljung');
box_test_res_1; 1-pchisq(box_test_res_1$statistic, df = 9); 1-pchisq(box_test_res_1$statistic, df = 12)
	Box-Ljung test

data:  ewrtn_ma_mod_1$residuals
X-squared = 17.603, df = 12, p-value = 0.1283
X-squared: 0.0400701658576007
X-squared: 0.12829019017333
In [393]:
box_test_res_2 = Box.test(ewrtn_ma_mod_2$residuals, lag = 12, type = 'Ljung');
box_test_res_2; 1-pchisq(box_test_res_2$statistic, df = 9); 1-pchisq(box_test_res_2$statistic, df = 12)
	Box-Ljung test

data:  ewrtn_ma_mod_2$residuals
X-squared = 17.621, df = 12, p-value = 0.1277
X-squared: 0.039839650376525
X-squared: 0.127706313281458
In [395]:
box_test_res_3 = Box.test(ewrtn_ma_mod_3$residuals, lag = 12, type = 'Ljung');
box_test_res_3; 1-pchisq(box_test_res_3$statistic, df = 9); 1-pchisq(box_test_res_3$statistic, df = 12)
	Box-Ljung test

data:  ewrtn_ma_mod_3$residuals
X-squared = 17.604, df = 12, p-value = 0.1283
X-squared: 0.0400560859710883
X-squared: 0.128254552132693

MA(9) forecasting (Sec. 2.5.4)

  • Plot the forecasting
In [503]:
eotr = 986
h = 10
npts = 12
freq = 12
order = c(0,0,9)
fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA)
seasonal = NULL
ewrtn_fc_res = plot_arima_forecast_fig(
    da_ts=ewrtn_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', transform.pars=F,
    main="Forecasts from ARIMA(0,0,9) with non-zero mean for ewrtn (Fig. 2.7)", 
    xlab="Year", ylab="Return", ylim=c(-0.25, 0.25)
)
2007.083 ; 2009
No description has been provided for this image
  • Forecasting results
In [504]:
ewrtn_multistep_ahead_forecast_tb = comb_forecast_res(
    forecast_obj=ewrtn_fc_res, 
    da_ts=ewrtn_ts, 
    eotr=eotr, 
    freq=freq
)
ewrtn_multistep_ahead_forecast_tb
Forecast method: ARIMA(0,0,9) with non-zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, transform.pars = transform_pars, fixed = fixed, 
    method = method)

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1844    0  -0.1206    0    0    0    0    0  0.1218     0.0128
s.e.  0.0295    0   0.0338    0    0    0    0    0  0.0312     0.0027

sigma^2 estimated as 0.005066:  log likelihood = 1206.44,  aic = -2402.88

Error measures:
                        ME       RMSE        MAE      MPE     MAPE      MASE
Training set -1.916748e-05 0.07117457 0.04750136 94.10055 214.5367 0.6922456
                   ACF1
Training set 0.02241295

Forecasts:
         Point Forecast       Lo 80      Hi 80      Lo 95     Hi 95
Mar 2008    0.004284404 -0.08692947 0.09549828 -0.1352152 0.1437840
Apr 2008    0.013560426 -0.07919123 0.10631208 -0.1282910 0.1554118
May 2008    0.015025565 -0.07772609 0.10777722 -0.1268258 0.1568770
Jun 2008    0.014453940 -0.07894734 0.10785522 -0.1283910 0.1572989
Jul 2008    0.012047087 -0.08135419 0.10544836 -0.1307978 0.1548920
Aug 2008    0.001806804 -0.09159447 0.09520808 -0.1410381 0.1446517
Sep 2008    0.012211891 -0.08118938 0.10561317 -0.1306330 0.1550568
Oct 2008    0.005515984 -0.08788529 0.09891726 -0.1373289 0.1483609
Nov 2008    0.008514008 -0.08488727 0.10191528 -0.1343309 0.1513589
Dec 2008    0.012792573 -0.08126721 0.10685236 -0.1310595 0.1566446
A Time Series: 1 × 10
MarAprMayJunJulAugSepOctNovDec
20080.0042844040.0135604260.0150255650.0144539400.0120470870.0018068040.0122118910.0055159840.0085140080.012792573
A Time Series: 1 × 10
MarAprMayJunJulAugSepOctNovDec
20080.071174570.072374500.072374500.072881400.072881400.072881400.072881400.072881400.072881400.07339524
A Time Series: 1 × 10
MarAprMayJunJulAugSepOctNovDec
2008-0.025971 0.031213 0.032178-0.087065-0.000994 0.014066-0.120888-0.206004-0.136597 0.043083
A Time Series: 10 × 3
ForecastStd. ErrorActual
Mar 20080.0042844040.07117457-0.025971
Apr 20080.0135604260.07237450 0.031213
May 20080.0150255650.07237450 0.032178
Jun 20080.0144539400.07288140-0.087065
Jul 20080.0120470870.07288140-0.000994
Aug 20080.0018068040.07288140 0.014066
Sep 20080.0122118910.07288140-0.120888
Oct 20080.0055159840.07288140-0.206004
Nov 20080.0085140080.07288140-0.136597
Dec 20080.0127925730.07339524 0.043083
In [ ]:

Find Equal-weighted monthly and yearly expected return

In [80]:
# µ, sigma, monthly pr ratio, yearly pr ratio
lewrtn_ts = ts(log(1+ewrtn_ts), frequency = 12, start = c(1926, 1))
ewrtn_basic_stats = apply(as.data.frame(lewrtn_ts), 2, basicStats)
ewrtn_basic_stats$x['Mean',];
ewrtn_basic_stats$x['Stdev',];
mnthly_ew_pr_ratio = exp(ewrtn_basic_stats$x['Mean',]+ewrtn_basic_stats$x['Stdev',]^2/2);
mnthly_ew_pr_ratio; mnthly_ew_pr_ratio^12;
0.009567
0.071387
1.01218873638879
1.15647966882156

Sec. 2.6

3M log price (Sec. 2.6.3)

In [82]:
da = read.table("../AFTS_data/Ch02/m-3m4608.txt", header = T)
da[1:5,]
A data.frame: 5 × 2
datertn
<int><dbl>
119460228-0.077922
219460330 0.018592
319460430-0.100000
419460531 0.209877
519460628 0.005128
In [83]:
lrtn = log(1+da$rtn)
lrtn_ts = ts(lrtn, frequency = 12, start = c(1946, 2))
In [516]:
plot_time_fig(lrtn_ts, main = "3M stock (Fig. 2.9)", xlab = "Year", ylab = "log(Return)")
No description has been provided for this image
In [517]:
par(bg = "white")
acf(lrtn_ts, main = "3M stock (Fig. 2.9)")
No description has been provided for this image
  • Print EACF results (Table 2.5)
In [522]:
find("eacf")
'package:TSA'
In [566]:
eacf_res <- eacf(lrtn_ts, ar.max = 6, ma.max = 12)
eacf_stats_tb = format(as.data.frame(eacf_res$eacf), digits = 3)
names(eacf_stats_tb) <- seq(from = 0, to = 12)
display(eacf_stats_tb)
display(eacf_res$symbol)
display(2/sqrt(length(lrtn_ts)))
AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12
0 o o x o o x o o o x o  x  o 
1 x o x o o x o o o o o  x  o 
2 x x x o o x o o o o o  o  o 
3 x x x o o o o o o o o  o  o 
4 x o x o o o o o o o o  o  o 
5 x x x o x o o o o o o  o  o 
6 x x x x x o o o o o o  o  o 
A data.frame: 7 × 13
0123456789101112
<I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>>
1-0.056-0.03796-0.0822-0.00465 0.017744 0.08207 0.008000.01267-0.030143-0.077770.04877 0.09091-0.0109
2-0.474 0.00958-0.0736-0.02090 0.001963 0.07722-0.028760.00259-0.006832-0.069370.03722 0.09383-0.0242
3-0.383-0.34759-0.0737 0.01595-0.005528 0.07720 0.026940.01197 0.000449-0.026780.02209 0.04284 0.0418
4-0.177 0.13811 0.3838-0.02240 0.002324 0.04192-0.023190.01542-0.004405-0.025370.01849 0.01002 0.0433
5 0.421 0.02873 0.4542-0.00787 0.000711 0.00254-0.014020.03052 0.011587 0.004220.01910-0.00434 0.0133
6-0.114 0.21353 0.4490 0.00958 0.202415-0.00627-0.003790.04032-0.012936-0.012330.03150 0.01172 0.0277
7-0.208-0.25037 0.2431 0.31109 0.167446-0.03875-0.003380.04292-0.010090-0.026000.00782 0.01056 0.0373
A matrix: 7 × 13 of type chr
0123456789101112
0ooxooxoooxoxo
1xoxooxoooooxo
2xxxooxooooooo
3xxxoooooooooo
4xoxoooooooooo
5xxxoxoooooooo
6xxxxxoooooooo
0.0727874525246839

Find 3M monthly and yearly expected return

  • For the log(Return)
    • $\mu_{monthly}\approx 0.010299$
    • $\sigma_{monthly}\approx 0.063719$
In [84]:
basic_3m_stats <- apply(as.data.frame(lrtn), 2, basicStats)
basic_3m_stats
# µ, std, and std of mean
basic_3m_stats$lrtn['Mean',]; basic_3m_stats$lrtn['Stdev',]; basic_3m_stats$lrtn['Stdev',]/sqrt(basic_3m_stats$lrtn['nobs',])
# 1 + Monthly simple grow rate
mnthly_3m_pr_ratio = exp(basic_3m_stats$lrtn['Mean',]+basic_3m_stats$lrtn['Stdev',]^2/2) 
mnthly_3m_pr_ratio; mnthly_3m_pr_ratio^12
$lrtn =
A data.frame: 16 × 1
X..newX..i
<dbl>
nobs755.000000
NAs 0.000000
Minimum -0.326128
Maximum 0.229484
1. Quartile -0.029969
3. Quartile 0.050031
Mean 0.010299
Median 0.008830
Sum 7.776056
SE Mean 0.002319
LCL Mean 0.005747
UCL Mean 0.014852
Variance 0.004060
Stdev 0.063719
Skewness -0.076466
Kurtosis 1.250917
0.010299
0.063719
0.00231897184371017
1.01240537159774
1.15945337518438

Sec. 2.7

3M log price

In [26]:
da = read.table("../AFTS_data/Ch02/m-3m4608.txt", header = T)
da[1:5,]
A data.frame: 5 × 2
datertn
<int><dbl>
119460228-0.077922
219460330 0.018592
319460430-0.100000
419460531 0.209877
519460628 0.005128
  • See the t-statistics and p-value of mean of log returns
In [27]:
lrtn = log(1+da$rtn)
lrtn_ts = ts(lrtn, frequency = 12, start = c(1946, 2))
In [28]:
t_stats = mean(lrtn)/(sd(lrtn)/sqrt(length(lrtn)))
pv = 2*(1-pnorm(abs(t_stats)))

mean(lrtn); 
sd(lrtn); 
sd(lrtn)/sqrt(length(lrtn)); # variance of the mean of the log returns
t_stats; 
pv
0.0102994120244493
0.063719098418256
0.00231897542551723
4.44136316026377
8.93907848342756e-06
  • Construct two price series
    • Using log return
    • Using mean-corrected log return
In [45]:
lgPr1 = cumsum(lrtn)
lgPr2 = cumsum(lrtn-mean(lrtn))
lgPr1_ts = ts(lgPr1, frequency = 12, start = c(1946, 1))
lgPr2_ts = ts(lgPr2, frequency = 12, start = c(1946, 1))
ln_ts = ts(mean(lrtn)*seq(from = 1, to = length(lrtn), by = 1), frequency = 12, start = c(1946, 1))
par(bg = "white")
plot(ln_ts, col = 'black', type = 'l', main = '3M constructed log-pr series (Fig. 2.10)', xlab = "Year", ylab = "log(Price)")
# plot(lgPr1_ts, type = 'p', pch = 1, main = '3M constructed log-pr series', xlab = "Year", ylab = "Log Price")
points(lgPr1_ts, col = 'red', pch = 1)
points(lgPr2_ts, col = 'blue', pch = 3)
# Add legend
legend(
    "topleft", 
    legend = c("y=0.0103*t", "sum(r_t)", "sum(r_t-mean(t_t))"), 
    col = c("black", "red", "blue"), 
    lty = c(1, NA, NA),
    pch = c(NA, 1, 3)
)
No description has been provided for this image
  • See some basic stats for 3M return
    • If we model the price series as a Geometric Brownian Motion, then the return series, i.e. diff(log(Price)), (with the constant time-step) should be sum of a constant mean and an iid Gaussian noise.
    • The Geometric Brownian Motion can be seen as an unit-root process in time-series.
In [46]:
pr_ts = ts(exp(lgPr1), frequency = 12, start = c(1946, 1))
plot_time_fig(ts = pr_ts, main = "Price", xlab = "Year", ylab = "Price/Price_0")
No description has been provided for this image
In [ ]:

Unit-root test (U.S. quarterly GDP)

In [3]:
da = read.table("../AFTS_data/Ch02/q-gdp4708.txt", header = T)
da[1:5,]
A data.frame: 5 × 4
yearmondaygdp
<int><int><int><dbl>
11947 11237.2
21947 41240.5
31947 71244.6
41947101254.4
51948 11260.4
In [7]:
lgdp = log(da$gdp)
m1 = ar(diff(lgdp), method = 'mle')
m1$order
10
In [570]:
adfTest(lgdp, lags = 10, type = c("c"))
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 10
  STATISTIC:
    Dickey-Fuller: -1.6109
  P VALUE:
    0.4569 

Description:
 Sat Feb 17 20:09:22 2024 by user: 
In [13]:
apply(as.data.frame(lgdp), 2, basicStats)
$lgdp =
A data.frame: 16 × 1
X..newX..i
<dbl>
nobs 248.000000
NAs 0.000000
Minimum 5.468904
Maximum 9.575872
1. Quartile 6.377224
3. Quartile 8.799009
Mean 7.612532
Median 7.664212
Sum1887.907986
SE Mean 0.081566
LCL Mean 7.451879
UCL Mean 7.773185
Variance 1.649936
Stdev 1.284498
Skewness -0.063967
Kurtosis -1.433830
In [600]:
lgdp_ts = ts(lgdp, frequency = 4, start = c(1947 ,1))
plot_unit_root_figs(
    lgdp_ts, freq = 4, start = c(1947, 2), lag_max = 20,
    main = 'U.S. quarterly GDP (Fig. 2.11)', xlab = "Year", ylab = "ln(GDP)"
)
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Full ARIMA workflow (S&P 500)

  • Unit-root testing
  • Order determination and parameter estimation
  • Box-Ljung test
  • Plot forcasting
In [207]:
# install.packages("RQuantLib")
also installing the dependencies ‘Rcpp’, ‘zoo’


Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

In [209]:
load_quantlib_calendars(ql_calendars = NULL, from = 1950, to = 2008, financial = TRUE)
Calendar QuantLib/TARGET loaded

Calendar QuantLib/Argentina loaded

Calendar QuantLib/Australia loaded

Calendar QuantLib/Brazil loaded

Calendar QuantLib/Canada loaded

Calendar QuantLib/Canada/Settlement loaded

Calendar QuantLib/Canada/TSX loaded

Calendar QuantLib/China loaded

Calendar QuantLib/CzechRepublic loaded

Calendar QuantLib/Denmark loaded

Calendar QuantLib/Finland loaded

Calendar QuantLib/Germany loaded

Calendar QuantLib/Germany/FrankfurtStockExchange loaded

Calendar QuantLib/Germany/Settlement loaded

Calendar QuantLib/Germany/Xetra loaded

Calendar QuantLib/Germany/Eurex loaded

Calendar QuantLib/HongKong loaded

Calendar QuantLib/Hungary loaded

Calendar QuantLib/Iceland loaded

Calendar QuantLib/India loaded

Calendar QuantLib/Indonesia loaded

Calendar QuantLib/Italy loaded

Calendar QuantLib/Italy/Settlement loaded

Calendar QuantLib/Italy/Exchange loaded

Calendar QuantLib/Japan loaded

Calendar QuantLib/Mexico loaded

Calendar QuantLib/NewZealand loaded

Calendar QuantLib/Norway loaded

Calendar QuantLib/Poland loaded

Calendar QuantLib/Russia loaded

Calendar QuantLib/SaudiArabia loaded

Calendar QuantLib/Singapore loaded

Calendar QuantLib/Slovakia loaded

Calendar QuantLib/SouthAfrica loaded

Calendar QuantLib/SouthKorea loaded

Calendar QuantLib/SouthKorea/KRX loaded

Calendar QuantLib/Sweden loaded

Calendar QuantLib/Switzerland loaded

Calendar QuantLib/Taiwan loaded

Calendar QuantLib/Turkey loaded

Calendar QuantLib/Ukraine loaded

Calendar QuantLib/UnitedKingdom loaded

Calendar QuantLib/UnitedKingdom/Settlement loaded

Calendar QuantLib/UnitedKingdom/Exchange loaded

Calendar QuantLib/UnitedKingdom/Metals loaded

Calendar QuantLib/UnitedStates loaded

Calendar QuantLib/UnitedStates/Settlement loaded

Calendar QuantLib/UnitedStates/NYSE loaded

Calendar QuantLib/UnitedStates/GovernmentBond loaded

Calendar QuantLib/UnitedStates/NERC loaded

In [15]:
da = read.table("../AFTS_data/Ch02/d-sp55008.txt", header = T)
da[1:10,]; dim(da)
A data.frame: 10 × 9
yearmondayopenhighlowclosevolumeadjclose
<int><int><int><dbl><dbl><dbl><dbl><dbl><dbl>
119501 316.6616.6616.6616.66126000016.66
219501 416.8516.8516.8516.85189000016.85
319501 516.9316.9316.9316.93255000016.93
419501 616.9816.9816.9816.98201000016.98
519501 917.0817.0817.0817.08252000017.08
6195011017.0317.0317.0317.03216000017.03
7195011117.0917.0917.0917.09263000017.09
8195011216.7616.7616.7616.76297000016.76
9195011316.6716.6716.6716.67333000016.67
10195011616.7216.7216.7216.72146000016.72
  1. 14662
  2. 9

Unit-root testing

In [16]:
lsp5 = log(da$close)
par(bg = "white")
lsp5_ts = ts(lsp5, frequency = 252, start = c(1950 , 1, 3))
plot(lsp5_ts, col = 'black', type = 'l', main = 'S&P 500 (Fig. 2.12)', xlab = "Year", ylab = "log(S&P 500)")
No description has been provided for this image
In [577]:
m2 = ar(diff(lsp5), method = 'mle')
m2$order
2
In [578]:
adfTest(lsp5, lags = 2, type = ("ct"))
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 2
  STATISTIC:
    Dickey-Fuller: -2.0179
  P VALUE:
    0.5708 

Description:
 Sat Feb 17 20:49:04 2024 by user: 
In [580]:
log(length(lsp5))
9.59301439178353
In [188]:
adfTest(lsp5, lags = 15, type = ("ct"))
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 15
  STATISTIC:
    Dickey-Fuller: -1.9946
  P VALUE:
    0.5807 

Description:
 Wed Feb 14 03:52:33 2024 by user: 

Order determination and parameter estimations

- With one-differencing
In [588]:
diff_lsp5_ts = ts(diff(lsp5), frequency = 252, start = c(1950 , 1, 4))
length(lsp5); length(diff_lsp5_ts);
display(head(lsp5, 10))
display(head(diff_lsp5_ts, 10))
14662
14661
  1. 2.8130106367387
  2. 2.82435065679837
  3. 2.82908719614504
  4. 2.8320361808832
  5. 2.83790818836042
  6. 2.8349764946746
  7. 2.8384934971275
  8. 2.81899509505394
  9. 2.8136106967627
  10. 2.81660560765656
A Time Series:
  1. 0.011340020059674
  2. 0.0047365393466734
  3. 0.00294898473815719
  4. 0.0058720074772225
  5. -0.0029316936858268
  6. 0.00351700245290232
  7. -0.0194984020735625
  8. -0.00538439829123405
  9. 0.00299491089385251
  10. 0.00833834491715413
In [601]:
plot_unit_root_figs(
    lsp5_ts, freq = 252, start = c(1950 , 1, 4), lag_max = 20,
    main = 'S&P500 Index', xlab = "Year", ylab = "ln(S&P500)"
)
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [602]:
diff_lsp5_ts = ts(diff(lsp5), frequency = 252, start = c(1950 , 1, 4))
In [605]:
par(bg = 'white')
pacf(diff_lsp5_ts)
No description has been provided for this image
In [606]:
par(bg = 'white')
acf(diff_lsp5_ts)
No description has been provided for this image
In [607]:
sp5_eacf_res <- eacf(diff_lsp5_ts, ar.max = 6, ma.max = 12)
sp5_eacf_stats_tb = format(as.data.frame(sp5_eacf_res$eacf), digits = 3)
names(sp5_eacf_stats_tb) <- seq(from = 0, to = 12)
display(sp5_eacf_stats_tb)
display(sp5_eacf_res$symbol)
# pp67, asymptotic standard error of EACF
display(2/sqrt(length(diff_lsp5_ts)))
AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12
0 x x o o o x o o o o o  x  o 
1 x x o o o o o o o o o  x  o 
2 x o o o o o o o o o o  x  o 
3 x x o o o o o o o o o  x  o 
4 x o o x o o o o o o o  x  o 
5 o x x x x x o o o o o  o  o 
6 x x x x x x o o o o o  o  o 
A data.frame: 7 × 13
0123456789101112
<I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>>
1 0.0694-0.0337-0.01280-0.002520 5.63e-04-0.01801-0.01383 0.000595-3.82e-03 0.002069-0.0064950.02531 0.000255
2 0.4254-0.0550-0.00613-0.004956-6.03e-05-0.01229-0.01443-0.001431-2.37e-03-0.001021-0.0001530.02540-0.000565
3-0.1919-0.0147 0.01549-0.000335 1.72e-03-0.01139-0.01638-0.000969 5.14e-03-0.000838-0.0002150.02057 0.006086
4-0.2794-0.1311-0.00067-0.007149 1.36e-02-0.00567-0.01635-0.000868 6.89e-04-0.001961-0.0003200.02061 0.004844
5 0.0655-0.0129 0.00773-0.194153-5.89e-03-0.00823-0.01319-0.002503 2.39e-05-0.001603-0.0006600.01698 0.009937
6 0.0085-0.1240-0.34471-0.230192 2.87e-02-0.04058-0.00963 0.003877-4.63e-03 0.004135-0.0081890.01272-0.007121
7-0.4457 0.2788-0.31801 0.262998-3.33e-01 0.11436 0.00744-0.008721-4.15e-03 0.002669-0.0012490.00202-0.001039
A matrix: 7 × 13 of type chr
0123456789101112
0xxoooxoooooxo
1xxoooooooooxo
2xooooooooooxo
3xxoooooooooxo
4xooxoooooooxo
5oxxxxxooooooo
6xxxxxxooooooo
0.0165176476943815
In [611]:
diff_sp5_mod = arima(diff_lsp5_ts, order = c(1,0,2), include.mean = F)
diff_sp5_mod
Call:
arima(x = diff_lsp5_ts, order = c(1, 0, 2), include.mean = F)

Coefficients:
         ar1     ma1      ma2
      0.0726  0.0001  -0.0373
s.e.  0.2893  0.2894   0.0251

sigma^2 estimated as 8.077e-05:  log likelihood = 48279.03,  aic = -96552.06
In [610]:
sp5_mod = arima(lsp5_ts, order = c(1,1,2), include.mean = T)
sp5_mod
Call:
arima(x = lsp5_ts, order = c(1, 1, 2), include.mean = T)

Coefficients:
         ar1     ma1      ma2
      0.0726  0.0001  -0.0373
s.e.  0.2893  0.2894   0.0251

sigma^2 estimated as 8.077e-05:  log likelihood = 48279.03,  aic = -96552.06

Box-Ljung test

  • pp33 very briefly talks about how to choose lag in the Box.test
  • Null hypothesis cannot be reject, so that the model is adequte.
In [615]:
diff_sp5_box_test = Box.test(diff_sp5_mod$residuals, lag = 5, type = 'Ljung')
diff_sp5_box_test
	Box-Ljung test

data:  diff_sp5_mod$residuals
X-squared = 1.6681, df = 5, p-value = 0.8929
In [616]:
diff_sp5_box_test_1 = Box.test(diff_sp5_mod$residuals, lag = 12, type = 'Ljung')
diff_sp5_box_test_1
	Box-Ljung test

data:  diff_sp5_mod$residuals
X-squared = 20.642, df = 12, p-value = 0.05587
In [612]:
sp5_box_test = Box.test(sp5_mod$residuals, lag = 5, type = 'Ljung')
sp5_box_test
	Box-Ljung test

data:  sp5_mod$residuals
X-squared = 1.6646, df = 5, p-value = 0.8933
In [613]:
sp5_box_test_1 = Box.test(sp5_mod$residuals, lag = 12, type = 'Ljung')
sp5_box_test_1
	Box-Ljung test

data:  sp5_mod$residuals
X-squared = 20.656, df = 12, p-value = 0.05565
In [ ]:

Plot forcasting

In [619]:
diff_sp5_mod$nobs
14661
In [676]:
npts = 22
eotr = 14661-npts
h = npts
freq = 252
order_1 = c(1,0,2)
fixed = NULL
seasonal = NULL
diff_sp5_fc_res = plot_arima_forecast_fig(
    da_ts=diff_lsp5_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order_1, seasonal=seasonal, fixed=fixed, method='ML', transform.pars=TRUE,
    main="Forecasts from ARIMA(1,0,2) with\n non-zero mean for ln(Return) of s&p500", 
    xlab="Year", ylab="ln(Return)", ylim=c(-0.05, 0.05)
)
[1] 14639
2008 ; 2008.179
No description has been provided for this image
In [683]:
order_2 = c(1,1,2)
sp5_fc_res = plot_arima_forecast_fig(
    da_ts=lsp5_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order_2, seasonal=seasonal, fixed=fixed, method='ML', transform.pars=TRUE,
    main="Forecasts from ARIMA(1,1,2) with\n non-zero mean for ln(Price) of s&p500", 
    xlab="Year", ylab="ln(Return)", ylim=c(7, 7.3)
)
[1] 14638
2008 ; 2008.183
No description has been provided for this image

Find S&P500 monthly and yearly expected return

  • We can think of
  • For the log(Return)
    • $\mu_{daily}\approx 0.000299$
    • $\sigma_{daily}\approx 0.009011$
In [85]:
# µ, std, and std of mean; daily, monthly, yearly pr-ratio
sp5_basic_stats <- apply(as.data.frame(diff(lsp5)), 2, basicStats)
display(sp5_basic_stats)
daily_exp_pr_ratio = exp(sp5_basic_stats$`diff(lsp5)`['Mean',]+sp5_basic_stats$`diff(lsp5)`['Stdev',]^2/2)
mnthly_exp_pr_ratio = daily_exp_pr_ratio^21
yrly_exp_pr_ratio = daily_exp_pr_ratio^252
sp5_basic_stats$`diff(lsp5)`['Mean',]; 
sp5_basic_stats$`diff(lsp5)`['Stdev',];
daily_exp_pr_ratio; mnthly_exp_pr_ratio; yrly_exp_pr_ratio;
$`diff(lsp5)` =
A data.frame: 16 × 1
X..newX..i
<dbl>
nobs14661.000000
NAs 0.000000
Minimum -0.228997
Maximum 0.087089
1. Quartile -0.004067
3. Quartile 0.004845
Mean 0.000299
Median 0.000443
Sum 4.382049
SE Mean 0.000074
LCL Mean 0.000153
UCL Mean 0.000445
Variance 0.000081
Stdev 0.009011
Skewness -1.276328
Kurtosis 33.590044
0.000299
0.009011
1.00033965673079
1.00715707054837
1.08934757636864

Sec. 2.8

J&J quarterly earnings

  • Starting at 1960 first quarter.
In [3]:
da = read.table("../AFTS_data/Ch02/q-jnj.txt", header = F) # needs to set `header = F`
head(da); length(da$V1)
A data.frame: 6 × 1
V1
<dbl>
10.71
20.63
30.85
40.44
50.61
60.69
84

Plot the acf of JnJ quarterly earnings

In [4]:
freq = 4
jnj_er = da$V1
jnj_er_ts = ts(jnj_er, frequency = freq, start = c(1960, 1))
jnj_lg_er = log(jnj_er)
jnj_lg_er_ts = ts(jnj_lg_er, frequency = freq, start = c(1960, 1))
jnj_lg_er_ts
A Time Series: 21 × 4
Qtr1Qtr2Qtr3Qtr4
1960-0.34249031-0.46203546-0.16251893-0.82098055
1961-0.49429632-0.37106368-0.08338161-0.59783700
1962-0.32850407-0.26136476-0.08338161-0.51082562
1963-0.18632958-0.22314355 0.00000000-0.26136476
1964-0.08338161 0.00000000 0.21511138 0.00000000
1965 0.14842001 0.26236426 0.37156356 0.22314355
1966 0.23111172 0.32208350 0.62057649 0.44468582
1967 0.42526774 0.46373402 0.60431597 0.62057649
1968 0.42526774 0.72754861 0.85015093 0.81093022
1969 0.77010822 0.88789126 0.99325177 0.81093022
1970 1.02604160 1.22964055 1.30562646 1.28093385
1971 1.28093385 1.46325540 1.46325540 1.39871688
1972 1.58103844 1.61740608 1.61740608 1.48387469
1973 1.71918878 1.76644166 1.88251383 1.66959184
1974 1.79674701 1.85473427 1.93585981 1.76644166
1975 1.93585981 2.04640169 2.05796251 1.81156210
1976 2.04640169 2.18717424 2.11384297 1.92278773
1977 2.25549349 2.32825284 2.25549349 2.16676526
1978 2.47485631 2.48989419 2.49732917 2.18717424
1979 2.64191040 2.56186769 2.69799987 2.30158459
1980 2.78501124 2.68580459 2.77383794 2.45186680
In [7]:
plot_time_fig(jnj_er_ts, main = "JnJ quarterly earnings (Fig. 2.13 (a))", xlab = "Year", ylab = "Earning")
plot_time_fig(jnj_lg_er_ts, main = "JnJ quarterly earnings (Fig. 2.13 (b))", xlab = "Year", ylab = "log(Earning)")
No description has been provided for this image
No description has been provided for this image
In [8]:
plot_pacf_acf(jnj_lg_er, main = "lg(earning) (Fig. 2.14 (a))")
plot_pacf_acf(diff(jnj_lg_er), main = "diff(lg(earning)) (Fig. 2.14 (b))")
plot_pacf_acf(diff(jnj_lg_er, lag = freq), main = "diff(lg(earning), 4) (Fig. 2.14 (c))")
plot_pacf_acf(diff(diff(jnj_lg_er), lag = freq), main = "diff(diff(lg(earning)), 4) (Fig. 2.14 (d))")
plot_pacf_acf(diff(diff(jnj_lg_er, lag = freq)), main = "diff(diff(lg(earning), 4)) (Fig. 2.14 (d)\')")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Apply airline model to the JnJ data.

In [7]:
find('arima')
  1. 'package:TSA'
  2. 'package:stats'
In [8]:
multi_seas_mod_1 <- arima(
    x = diff(diff(jnj_lg_er_ts), lag = freq),
    order = c(0,0,1),
    seasonal = list(order = c(0,0,1), period = freq),
    method = 'ML',
    include.mean = F
)
summary(multi_seas_mod_1); sqrt(multi_seas_mod_1$sigma2)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = diff(diff(jnj_lg_er_ts), lag = freq), order = c(0, 0, 1), seasonal = list(order = c(0, 
    0, 1), period = freq), include.mean = F, method = "ML")

Coefficients:
          ma1     sma1
      -0.6809  -0.3146
s.e.   0.0982   0.1070

sigma^2 estimated as 0.007931:  log likelihood = 78.38,  aic = -152.75

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.0890540593592921
In [9]:
multi_seas_mod_2 <- arima(
    x = diff(diff(jnj_lg_er_ts), lag = freq),
    order = c(0,0,1),
    seasonal = list(order = c(0,0,1), period = freq),
    method = 'CSS-ML',
    include.mean = F
)
summary(multi_seas_mod_2); sqrt(multi_seas_mod_2$sigma2)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = diff(diff(jnj_lg_er_ts), lag = freq), order = c(0, 0, 1), seasonal = list(order = c(0, 
    0, 1), period = freq), include.mean = F, method = "CSS-ML")

Coefficients:
          ma1     sma1
      -0.6809  -0.3146
s.e.   0.0982   0.1070

sigma^2 estimated as 0.007931:  log likelihood = 78.38,  aic = -152.75

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.0890540586073813
In [10]:
multi_seas_mod_3 <- arima(
    x = diff(diff(jnj_lg_er_ts), lag = freq),
    order = c(0,0,1),
    seasonal = list(order = c(0,0,1), period = freq),
    method = 'CSS',
    include.mean = F
)
summary(multi_seas_mod_3); sqrt(multi_seas_mod_3$sigma2)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = diff(diff(jnj_lg_er_ts), lag = freq), order = c(0, 0, 1), seasonal = list(order = c(0, 
    0, 1), period = freq), include.mean = F, method = "CSS")

Coefficients:
          ma1     sma1
      -0.5228  -0.2752
s.e.   0.0956   0.0992

sigma^2 estimated as 0.008665:  part log likelihood = 75.47

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.0930862384199885
In [123]:
# multi_seas_mod_4 <- arima(
#     x = diff(diff(jnj_lg_er_ts), lag = freq),
#     order = c(1,0,1),
#     seasonal = list(order = c(1,0,1), period = freq),
#     method = 'ML',
#     include.mean = F,
#     transform.pars = F
# )
# # multi_seas_mod_4 <- arima(
# #     x = diff(diff(jnj_lg_er_ts), lag = freq),
# #     order = c(1,0,1),
# #     seasonal = list(order = c(1,0,1), period = freq),
# #     fixed = c(1, NA, 1, NA), 
# #     method = 'ML',
# #     include.mean = F,
# #     transform.pars = F
# # )
# summary(multi_seas_mod_4); sqrt(multi_seas_mod_4$sigma2)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = diff(diff(jnj_lg_er_ts), lag = freq), order = c(1, 0, 1), seasonal = list(order = c(1, 
    0, 1), period = freq), include.mean = F, transform.pars = F, method = "ML")

Coefficients:
          ar1      ma1     sar1     sma1
      -0.0180  -0.6626  -0.2075  -0.1382
s.e.   0.2216   0.1815   0.3080   0.2950

sigma^2 estimated as 0.007889:  log likelihood = 78.57,  aic = -149.14

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.0888177229169977

Box-Ljung test

  • Cannot reject the null-hypothesis. So the model appears to be adequate.
In [121]:
box_test <- Box.test(multi_seas_mod_1$residuals, lag = 12, type = 'Ljung')
box_test
	Box-Ljung test

data:  multi_seas_mod_1$residuals
X-squared = 10.15, df = 12, p-value = 0.6028

Forecast

  • (0):

$$ \begin{eqnarray} z_t &=& (1-B)(1-B^4)x_t = (1-\theta B)(1-\beta B^4)a_t = a_t-\theta a_{t-1}-\beta a_{t-4}+\theta \beta a_{t-5} \\ y_t &=& (1-B^4)x_t \\ z_t &=& (1-B)y_t \end{eqnarray} $$

$$ \begin{eqnarray} y_{t+1}&=&z_{t+1}+y_t \\ x_{t+4}&=&y_{t+4}+x_t \end{eqnarray} $$

or

$$ \begin{eqnarray} y_{t}&=&z_{t}+y_{t-1} \\ x_{t}&=&y_{t}+x_{t-4} \end{eqnarray} $$

The expectation is

$$ \begin{eqnarray} E[y_{t+1}]&=&E[z_{t+1}]+E[y_t] \\ E[x_{t+4}]&=&E[y_{t+4}]+E[x_t] \end{eqnarray} $$

  • (1) (following are stationary results, but for forecasting needs to consider the time-zero):

$$ \begin{eqnarray} Var[z_t]&=&(1+\theta^2+\beta^2+\theta^2\beta^2)\sigma^2 \\ Cov(z_{t+1},z_t)&=&-\theta(1+\beta^2)\sigma^2 \\ Cov(z_{t+2},z_t)&=&0 \\ Cov(z_{t+3},z_t)&=&\theta\beta\sigma^2 \\ Cov(z_{t+4},z_t)&=&-\beta(1+\theta^2)\sigma^2 \\ Cov(z_{t+5},z_t)&=&\theta\beta\sigma^2 \\ Cov(z_{t+i},z_t)&=&0, \forall i>5 \end{eqnarray} $$

  • (2):

$$ \begin{eqnarray} Var[y_t]&=&Var[z_t]+Var[y_{t-1}]+Cov(z_t, y_{t-1}) \\ Cov(y_{t+1}, y_t)&=&Cov(z_{t+1}, y_t)+Var[y_t]=Cov(z_{t+1},z_t+z_{t-1}+z_{t-2}+z_{t-3}+z_{t-4})+Var[y_t] \\ Cov(y_{t+2}, y_t)&=&Cov(z_{t+2}, y_t)+Cov(y_{t+1}, y_t)=Cov(z_{t+2},z_t+z_{t-1}+z_{t-2}+z_{t-3})+Cov(y_{t+1}, y_t) \\ Cov(y_{t+3}, y_t)&=&Cov(z_{t+3}, y_t)+Cov(y_{t+2}, y_t)=Cov(z_{t+3},z_t+z_{t-1}+z_{t-2})+Cov(y_{t+2}, y_t) \\ Cov(y_{t+4}, y_t)&=&Cov(z_{t+4}, y_t)+Cov(y_{t+3}, y_t)=Cov(z_{t+4},z_t+z_{t-1})+Cov(y_{t+3}, y_t) \\ Cov(y_{t+5}, y_t)&=&Cov(z_{t+5}, y_t)+Cov(y_{t+4}, y_t)=Cov(z_{t+5},z_t)+Cov(y_{t+4}, y_t) \\ Cov(y_{t+i}, y_t)&=&Cov(y_{t+5}, y_t), \forall i>5 \end{eqnarray} $$

  • (3):

$$ \begin{eqnarray} Var[x_t]&=&Var[y_t]+Var[x_{t-4}]+Cov(y_t,y_{t-4}+y_{t-8}+\cdots) \end{eqnarray} $$

In [227]:
# plot_arima_forecast_fig <- function(
#     da_ts, eotr, h, npts, frequency,
#     order, seasonal, fixed, method, include.mean, transform.pars,
#     main, xlab, ylab, ylim = NULL) {
#     par(bg = "white")
#     # arima model
#     tr_da_ts <- ts(da_ts[1:eotr], frequency = frequency, start = start(da_ts))
#     if (is.null(transform.pars)) {
#         ts_fm3 <- arima(tr_da_ts, order = order, fixed = fixed, seasonal = seasonal, method = method, include.mean = include.mean)
#     } else {
#         ts_fm3 <- arima(tr_da_ts, order = order, fixed = fixed, seasonal = seasonal, method = method, include.mean = include.mean, transform.pars = transform.pars)
#     }
#     print(ts_fm3$nobs)
#     # Forecast
#     ts_fm3$x <- tr_da_ts # https://stackoverflow.com/a/42464130/4307919
#     ts_fc_res <- forecast(ts_fm3, h = h)
#     # Plot forecast
#     if (is.null(npts)) {
#         npts <- eotr
#     }
#     xmin <- time(da_ts)[eotr] - npts / frequency
#     xmax <- time(da_ts)[eotr] + (max(h, length(da_ts) - eotr) + 1) / frequency
#     cat(xmin, ";", xmax)

#     plot(ts_fc_res, xlim = c(xmin, xmax), ylim = ylim, main = main, xlab = xlab, ylab = ylab)
#     # Plot forecast mean (prepend the last observed data in the training dataset)
#     dummy_1st_fmean_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$mean)), frequency = frequency, start = end(tr_da_ts))
#     lines(dummy_1st_fmean_ts)
#     points(dummy_1st_fmean_ts, pch = 1)
#     # Plot confidence interval (95%)
#     dummy_1st_flower_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$lower[, 2])), frequency = frequency, start = end(tr_da_ts))
#     dummy_1st_fupper_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$upper[, 2])), frequency = frequency, start = end(tr_da_ts))
#     lines(dummy_1st_flower_ts, lty = 2)
#     lines(dummy_1st_fupper_ts, lty = 2)
#     # Plot original data
#     orig_plot_ts <- ts(da_ts[(eotr - npts + 1):length(da_ts)], frequency = frequency, start = time(da_ts)[eotr] - (npts - 1) / frequency)
#     lines(orig_plot_ts)
#     points(orig_plot_ts, pch = 19)
#     ts_fc_res
# }

# comb_forecast_res <- function(forecast_obj, da_ts, eotr, freq) {
#     display(summary(forecast_obj))
#     fc_std_err = (forecast_obj$upper[,2]-forecast_obj$lower[,2])/2/qnorm(p = 0.975)
#     actual_ts = ts(da_ts[(eotr+1):length(da_ts)], frequency = freq, start = time(da_ts)[eotr+1])
#     display(forecast_obj$mean); display(fc_std_err); display(actual_ts)
#     multistep_ahead_forecast_tb = cbind(forecast_obj$mean, fc_std_err, actual_ts)
#     dimnames(multistep_ahead_forecast_tb)[[2]] <- c("Forecast", "Std. Error", "Actual")
#     multistep_ahead_forecast_tb
# }
In [229]:
npts = 6
h = 8
eotr = 79-h
freq = 4
order = c(0,0,1)
seasonal = list(order = order, period = freq)
fixed = NULL
include.mean = F
diff_seasonal_reg_jnj_ts = ts(diff(diff(jnj_lg_er_ts), lag = freq), frequency = freq, start = c(1961, 2))
diff_seasonal_reg_jnj_fc_res = plot_arima_forecast_fig(
    da_ts=diff_seasonal_reg_jnj_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', include.mean=include.mean, transform.pars=TRUE,
    main="Forecasts from ARIMA(0,0,1)(0,0,1)[4] with\n non-zero mean for (1-B)(1-B^4)ln(QtlyEarning) of JnJ", 
    xlab="Year", ylab="(1-B)(1-B^4)ln(QtlyEarning)"#, ylim=c(-0.05, 0.05)
)
summary(diff_seasonal_reg_jnj_fc_res)
[1] 71
1977.25 ; 1981
Forecast method: ARIMA(0,0,1)(0,0,1)[4] with zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.6559  -0.3492
s.e.   0.1094   0.1104

sigma^2 estimated as 0.008409:  log likelihood = 68.28,  aic = -132.57

Error measures:
                     ME      RMSE        MAE  MPE MAPE      MASE       ACF1
Training set 0.01083226 0.0916983 0.07697376 -Inf  Inf 0.5682199 0.05218151

Forecasts:
        Point Forecast       Lo 80     Hi 80       Lo 95     Hi 95
1979 Q1     0.09037218 -0.02714392 0.2078883 -0.08935319 0.2700975
1979 Q2     0.02639844 -0.11413996 0.1669369 -0.18853650 0.2413334
1979 Q3    -0.02343313 -0.16397154 0.1171053 -0.23836807 0.1915018
1979 Q4     0.06359381 -0.07694460 0.2041322 -0.15134113 0.2785288
1980 Q1    -0.03378421 -0.18019118 0.1126228 -0.25769434 0.1901259
1980 Q2     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
1980 Q3     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
1980 Q4     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
No description has been provided for this image
In [251]:
jnj_multistep_ahead_forecast_tb <- comb_forecast_res(diff_seasonal_reg_jnj_fc_res, diff_seasonal_reg_jnj_ts, eotr, freq)
jnj_multistep_ahead_forecast_tb
Forecast method: ARIMA(0,0,1)(0,0,1)[4] with zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.6559  -0.3492
s.e.   0.1094   0.1104

sigma^2 estimated as 0.008409:  log likelihood = 68.28,  aic = -132.57

Error measures:
                     ME      RMSE        MAE  MPE MAPE      MASE       ACF1
Training set 0.01083226 0.0916983 0.07697376 -Inf  Inf 0.5682199 0.05218151

Forecasts:
        Point Forecast       Lo 80     Hi 80       Lo 95     Hi 95
1979 Q1     0.09037218 -0.02714392 0.2078883 -0.08935319 0.2700975
1979 Q2     0.02639844 -0.11413996 0.1669369 -0.18853650 0.2413334
1979 Q3    -0.02343313 -0.16397154 0.1171053 -0.23836807 0.1915018
1979 Q4     0.06359381 -0.07694460 0.2041322 -0.15134113 0.2785288
1980 Q1    -0.03378421 -0.18019118 0.1126228 -0.25769434 0.1901259
1980 Q2     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
1980 Q3     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
1980 Q4     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
A Time Series: 2 × 4
Qtr1Qtr2Qtr3Qtr4
1979 0.09037218 0.02639844-0.02343313 0.06359381
1980-0.03378421 0.00000000 0.00000000 0.00000000
A Time Series: 2 × 4
Qtr1Qtr2Qtr3Qtr4
19790.09169830.10966270.10966270.1096627
19800.11424200.11615640.11615640.1161564
A Time Series: 2 × 4
Qtr1Qtr2Qtr3Qtr4
1979 0.14664510-0.09508059 0.12869720-0.08626034
1980 0.02869049-0.01916394-0.04809882 0.07444413
A Time Series: 8 × 3
ForecastStd. ErrorActual
1979 Q1 0.090372180.0916983 0.14664510
1979 Q2 0.026398440.1096627-0.09508059
1979 Q3-0.023433130.1096627 0.12869720
1979 Q4 0.063593810.1096627-0.08626034
1980 Q1-0.033784210.1142420 0.02869049
1980 Q2 0.000000000.1161564-0.01916394
1980 Q3 0.000000000.1161564-0.04809882
1980 Q4 0.000000000.1161564 0.07444413
  • Use DP to calculate forecasting expectation and variance. $$ \begin{eqnarray} E[y_{t+1}]&=&E[z_{t+1}]+E[y_t] \\ E[x_{t+4}]&=&E[y_{t+4}]+E[x_t] \end{eqnarray} $$
In [360]:
# Assume freq=s>=4
# z_t
diff_seasonal_reg_jnj_ts = ts(diff(diff(jnj_lg_er_ts), lag = freq), frequency = freq, start = c(1961, 2))
# y_t
diff_seasonal_jnj_ts = ts(diff(jnj_lg_er_ts, lag = freq), frequency = freq, start = c(1961, 1))
fstart = c(1979, 1)
iz = 71
iy = 72
ix = 76
h = 8
# E[y_t]
e_y <- list()
for (i in 1:h) {
    if (i==1) {
        e_y[[paste(i)]] = diff_seasonal_reg_jnj_fc_res$mean[i] + diff_seasonal_jnj_ts[length(diff_seasonal_jnj_ts)-h]
    } else {
        e_y[[paste(i)]] = diff_seasonal_reg_jnj_fc_res$mean[i] + e_y[[paste(i-1)]]
    }
}
# E[x_t]
e_x <- list()
for (i in 1:h) {
    if (i<=freq) {
        e_x[[paste(i)]] = e_y[[paste(i)]] + jnj_lg_er_ts[length(jnj_lg_er_ts)-h+i-freq]
    } else {
        e_x[[paste(i)]] = e_y[[paste(i)]] + e_x[[paste(i-freq)]]
    }
}
theta = -as.numeric(diff_seasonal_reg_jnj_fc_res$model$coef[1])
beta = -as.numeric(diff_seasonal_reg_jnj_fc_res$model$coef[2])
sigma2 = as.numeric(diff_seasonal_reg_jnj_fc_res$model$sigma2)
gen_key <- function(i,j) {paste("(",i,",",j,")", sep = "")}
# V[z_t]
# The following are unconditional, we need to calculate conditional
# v_z <- list(
#     "0"=as.numeric((1+theta^2+beta^2+theta^2*beta^2)*sigma2), 
#     "1"=as.numeric(-theta*(1+beta^2)*sigma2),
#     "2"=0,
#     "3"=as.numeric(theta*beta*sigma2),
#     "4"=as.numeric(-beta*(1+theta^2)*sigma2),
#     "5"=as.numeric(theta*beta*sigma2)
# )
v_z <- list()
for (i in 1:h) {
    for (j in 0:h) {
        v_z[[gen_key(i,j)]] = 0
        if (j==0) {
            v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] + sigma2
            if (i-1>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] + theta^2*sigma2}
            if (i-freq>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] + beta^2*sigma2}
            if (i-1-freq>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] + theta^2*beta^2*sigma2}
        } else if (j==1) {
            v_z[[gen_key(i,j)]] = -theta*sigma2
            if (i-freq>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] - theta*beta^2*sigma2}
        } else if (j==freq-1) {
            if (i-1>0) {v_z[[gen_key(i,j)]] = theta*beta*sigma2}
        } else if (j==freq) {
            v_z[[gen_key(i,j)]] = -beta*sigma2
            if (i-1>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] - beta*theta^2*sigma2}
        } else if (j==freq+1) {v_z[[gen_key(i,j)]] = theta*beta*sigma2}
        else {v_z[[gen_key(i,j)]] = 0}
    }
    for (j in 1:h) {
        v_z[[gen_key(i,-j)]] = 0
        if (j==1) {
            if (i-1>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] - theta*sigma2}
            if (i-1-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] - theta*beta^2*sigma2}
        } else if (j==freq-1) {
            if (i-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] + theta*beta*sigma2}
        } else if (j==freq) {
            if (i-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] - beta*sigma2}
            if (i-1-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] - beta*theta^2*sigma2}
        } else if (j==freq+1) {
            if (i-1-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] + theta*beta*sigma2}
        }
    }
}
# V[y_t]
v_tau <- list(
    "0"=sigma2,
    "1"=-theta*sigma2,
    "4"=-beta*sigma2,
    "5"=theta*beta*sigma2
)
v_y <- list()
for (i in 1:h) {
    v_y[[gen_key(i,0)]] = v_z[[gen_key(i,0)]]
    if (i-1>0) {v_y[[gen_key(i,0)]] = v_y[[gen_key(i,0)]] + v_y[[gen_key(i-1,0)]] - theta*v_tau[["0"]]}
    if (i-freq>0) {v_y[[gen_key(i,0)]] = v_y[[gen_key(i,0)]] - beta*(v_tau[["0"]]+v_tau[["1"]])}
    if (i-1-freq>0) {v_y[[gen_key(i,0)]] = v_y[[gen_key(i,0)]] + theta*beta*(v_tau[["0"]]+v_tau[["1"]]+v_tau[["4"]])}
    for (j in 1:h) {
        v_y[[gen_key(i,j)]] = v_y[[gen_key(i,j-1)]]
        if (j<=freq+1) {
            for (k in j:(freq+1)) {
                if (i+j-k>0) {v_y[[gen_key(i,j)]] = v_y[[gen_key(i,j)]] + v_z[[gen_key(i+j-k,k)]]}
                else {v_y[[gen_key(i,j)]] = v_y[[gen_key(i,j)]] + v_z[[gen_key(i+j,-k)]]}
            }
        }
    }
}
# V[x_t]
v_x <- list()
for (i in 1:h) {
    v_x[[paste(i)]] = v_y[[gen_key(i,0)]]
    if (i>freq) {v_x[[paste(i)]] = v_x[[paste(i)]] + v_x[[paste(i-freq)]]}
    j = freq
    while (TRUE) {
        if (i>j) {
            v_x[[paste(i)]] = v_x[[paste(i)]] + v_y[[gen_key(i-j,j)]]
            j = j + freq
        }
        else {break}
    }
}
In [361]:
e_y; e_x
$`1`
0.110781166566631
$`2`
0.137179611395488
$`3`
0.113746480815894
$`4`
0.177340290966718
$`5`
0.143556085705893
$`6`
0.143556085705893
$`7`
0.143556085705893
$`8`
0.143556085705893
$`1`
2.58563748050113
$`2`
2.62707380269453
$`3`
2.61107565060245
$`4`
2.36451453244944
$`5`
2.72919356620702
$`6`
2.77062988840042
$`7`
2.75463173630834
$`8`
2.50807061815533
In [416]:
# My calculation of std.error match the results from the package
for (i in 1:h) {
    cat(i, ":", sqrt(v_z[[gen_key(i,0)]]), "\n")
}
jnj_multistep_ahead_forecast_tb
1 : 0.0916983 
2 : 0.1096627 
3 : 0.1096627 
4 : 0.1096627 
5 : 0.114242 
6 : 0.1161564 
7 : 0.1161564 
8 : 0.1161564 
A Time Series: 8 × 3
ForecastStd. ErrorActual
1979 Q1 0.090372180.0916983 0.14664510
1979 Q2 0.026398440.1096627-0.09508059
1979 Q3-0.023433130.1096627 0.12869720
1979 Q4 0.063593810.1096627-0.08626034
1980 Q1-0.033784210.1142420 0.02869049
1980 Q2 0.000000000.1161564-0.01916394
1980 Q3 0.000000000.1161564-0.04809882
1980 Q4 0.000000000.1161564 0.07444413
In [363]:
v_z
$`(1,0)`
0.00840857826352835
$`(1,1)`
-0.00551512456267075
$`(1,2)`
0
$`(1,3)`
0
$`(1,4)`
-0.00293623482570831
$`(1,5)`
0.001925854799886
$`(1,6)`
0
$`(1,7)`
0
$`(1,8)`
0
$`(1,-1)`
0
$`(1,-2)`
0
$`(1,-3)`
0
$`(1,-4)`
0
$`(1,-5)`
0
$`(1,-6)`
0
$`(1,-7)`
0
$`(1,-8)`
0
$`(2,0)`
0.0120259078510645
$`(2,1)`
-0.00551512456267075
$`(2,2)`
0
$`(2,3)`
0.001925854799886
$`(2,4)`
-0.00419938880704877
$`(2,5)`
0.001925854799886
$`(2,6)`
0
$`(2,7)`
0
$`(2,8)`
0
$`(2,-1)`
-0.00551512456267075
$`(2,-2)`
0
$`(2,-3)`
0
$`(2,-4)`
0
$`(2,-5)`
0
$`(2,-6)`
0
$`(2,-7)`
0
$`(2,-8)`
0
$`(3,0)`
0.0120259078510645
$`(3,1)`
-0.00551512456267075
$`(3,2)`
0
$`(3,3)`
0.001925854799886
$`(3,4)`
-0.00419938880704877
$`(3,5)`
0.001925854799886
$`(3,6)`
0
$`(3,7)`
0
$`(3,8)`
0
$`(3,-1)`
-0.00551512456267075
$`(3,-2)`
0
$`(3,-3)`
0
$`(3,-4)`
0
$`(3,-5)`
0
$`(3,-6)`
0
$`(3,-7)`
0
$`(3,-8)`
0
$`(4,0)`
0.0120259078510645
$`(4,1)`
-0.00551512456267075
$`(4,2)`
0
$`(4,3)`
0.001925854799886
$`(4,4)`
-0.00419938880704877
$`(4,5)`
0.001925854799886
$`(4,6)`
0
$`(4,7)`
0
$`(4,8)`
0
$`(4,-1)`
-0.00551512456267075
$`(4,-2)`
0
$`(4,-3)`
0
$`(4,-4)`
0
$`(4,-5)`
0
$`(4,-6)`
0
$`(4,-7)`
0
$`(4,-8)`
0
$`(5,0)`
0.0130512268385915
$`(5,1)`
-0.00618762373618857
$`(5,2)`
0
$`(5,3)`
0.001925854799886
$`(5,4)`
-0.00419938880704877
$`(5,5)`
0.001925854799886
$`(5,6)`
0
$`(5,7)`
0
$`(5,8)`
0
$`(5,-1)`
-0.00551512456267075
$`(5,-2)`
0
$`(5,-3)`
0.001925854799886
$`(5,-4)`
-0.00293623482570831
$`(5,-5)`
0
$`(5,-6)`
0
$`(5,-7)`
0
$`(5,-8)`
0
$`(6,0)`
0.0134923140942493
$`(6,1)`
-0.00618762373618857
$`(6,2)`
0
$`(6,3)`
0.001925854799886
$`(6,4)`
-0.00419938880704877
$`(6,5)`
0.001925854799886
$`(6,6)`
0
$`(6,7)`
0
$`(6,8)`
0
$`(6,-1)`
-0.00618762373618857
$`(6,-2)`
0
$`(6,-3)`
0.001925854799886
$`(6,-4)`
-0.00419938880704877
$`(6,-5)`
0.001925854799886
$`(6,-6)`
0
$`(6,-7)`
0
$`(6,-8)`
0
$`(7,0)`
0.0134923140942493
$`(7,1)`
-0.00618762373618857
$`(7,2)`
0
$`(7,3)`
0.001925854799886
$`(7,4)`
-0.00419938880704877
$`(7,5)`
0.001925854799886
$`(7,6)`
0
$`(7,7)`
0
$`(7,8)`
0
$`(7,-1)`
-0.00618762373618857
$`(7,-2)`
0
$`(7,-3)`
0.001925854799886
$`(7,-4)`
-0.00419938880704877
$`(7,-5)`
0.001925854799886
$`(7,-6)`
0
$`(7,-7)`
0
$`(7,-8)`
0
$`(8,0)`
0.0134923140942493
$`(8,1)`
-0.00618762373618857
$`(8,2)`
0
$`(8,3)`
0.001925854799886
$`(8,4)`
-0.00419938880704877
$`(8,5)`
0.001925854799886
$`(8,6)`
0
$`(8,7)`
0
$`(8,8)`
0
$`(8,-1)`
-0.00618762373618857
$`(8,-2)`
0
$`(8,-3)`
0.001925854799886
$`(8,-4)`
-0.00419938880704877
$`(8,-5)`
0.001925854799886
$`(8,-6)`
0
$`(8,-7)`
0
$`(8,-8)`
0
In [364]:
v_y
$`(1,0)`
0.00840857826352835
$`(1,1)`
0.0028934537008576
$`(1,2)`
0.0028934537008576
$`(1,3)`
0.0028934537008576
$`(1,4)`
-4.27811248507096e-05
$`(1,5)`
0.00188307367503529
$`(1,6)`
0.00188307367503529
$`(1,7)`
0.00188307367503529
$`(1,8)`
0.00188307367503529
$`(2,0)`
0.0149193615519221
$`(2,1)`
0.00940423698925133
$`(2,2)`
0.00940423698925133
$`(2,3)`
0.00839385696342902
$`(2,4)`
0.00612032295626625
$`(2,5)`
0.00804617775615225
$`(2,6)`
0.00804617775615225
$`(2,7)`
0.00804617775615225
$`(2,8)`
0.00804617775615225
$`(3,0)`
0.0214301448403158
$`(3,1)`
0.015915020277645
$`(3,2)`
0.0149046402518227
$`(3,3)`
0.014556961044546
$`(3,4)`
0.0122834270373832
$`(3,5)`
0.0142092818372692
$`(3,6)`
0.0142092818372692
$`(3,7)`
0.0142092818372692
$`(3,8)`
0.0142092818372692
$`(4,0)`
0.0279409281287095
$`(4,1)`
0.0214154235402165
$`(4,2)`
0.0210677443329397
$`(4,3)`
0.0207200651256629
$`(4,4)`
0.0184465311185002
$`(4,5)`
0.0203723859183862
$`(4,6)`
0.0203723859183862
$`(4,7)`
0.0203723859183862
$`(4,8)`
0.0203723859183862
$`(5,0)`
0.0344666503788079
$`(5,1)`
0.0279313474353426
$`(5,2)`
0.0275836682280658
$`(5,3)`
0.027235989020789
$`(5,4)`
0.0249624550136263
$`(5,5)`
0.0268883098135123
$`(5,6)`
0.0268883098135123
$`(5,7)`
0.0268883098135123
$`(5,8)`
0.0268883098135123
$`(6,0)`
0.0414236615295919
$`(6,1)`
0.0348883585861266
$`(6,2)`
0.0345406793788498
$`(6,3)`
0.034193000171573
$`(6,4)`
0.0319194661644103
$`(6,5)`
0.0338453209642963
$`(6,6)`
0.0338453209642963
$`(6,7)`
0.0338453209642963
$`(6,8)`
0.0338453209642963
$`(7,0)`
0.0483806726803759
$`(7,1)`
0.0418453697369105
$`(7,2)`
0.0414976905296338
$`(7,3)`
0.041150011322357
$`(7,4)`
0.0388764773151942
$`(7,5)`
0.0408023321150802
$`(7,6)`
0.0408023321150802
$`(7,7)`
0.0408023321150802
$`(7,8)`
0.0408023321150802
$`(8,0)`
0.0553376838311598
$`(8,1)`
0.0488023808876945
$`(8,2)`
0.0484547016804177
$`(8,3)`
0.048107022473141
$`(8,4)`
0.0458334884659782
$`(8,5)`
0.0477593432658642
$`(8,6)`
0.0477593432658642
$`(8,7)`
0.0477593432658642
$`(8,8)`
0.0477593432658642
In [365]:
v_x
$`1`
0.00840857826352835
$`2`
0.0149193615519221
$`3`
0.0214301448403158
$`4`
0.0279409281287095
$`5`
0.0428324475174856
$`6`
0.0624633460377802
$`7`
0.0820942445580749
$`8`
0.101725143078369
In [366]:
theta; beta; sigma2
0.655892636046718
0.34919515923923
0.00840857826352835
  • Plot the forecasting
    • The resulting plot is somewhat different from the textbook Fig. 2.15.
In [384]:
e_qer_no_var_arr <- c()
e_lg_qer_arr <- c()
e_qer_arr <- c()
v_qer_arr <- c()
for (i in 1:h) {
    e_qer_no_var_arr <- c(e_qer_no_var_arr, exp(e_x[[paste(i)]]))
    e_lg_qer_arr <- c(e_lg_qer_arr, e_x[[paste(i)]])
    e_qer_arr <- c(e_qer_arr, exp(e_x[[paste(i)]]+v_x[[paste(i)]]/2))
    v_qer_arr <- c(v_qer_arr, v_x[[paste(i)]])
    cat(i, ":", 
        exp(e_x[[paste(i)]]), exp(e_x[[paste(i)]]+v_x[[paste(i)]]/2), 
        e_x[[paste(i)]], v_x[[paste(i)]], "\n",
        sep = ","
    )
}
1,:,13.27175,13.32766,2.585637,0.008408578,
2,:,13.83323,13.93681,2.627074,0.01491936,
3,:,13.61369,13.76034,2.611076,0.02143014,
4,:,10.63887,10.78855,2.364515,0.02794093,
5,:,15.32053,15.65217,2.729194,0.04283245,
6,:,15.96869,16.47529,2.77063,0.06246335,
7,:,15.71525,16.37374,2.754632,0.08209424,
8,:,12.28121,12.92202,2.508071,0.1017251,
In [3]:
help(legend)
In [415]:
par(bg = "white")
da_ts = jnj_er_ts
lg_da_ts = jnj_lg_er_ts
eotr = length(da_ts)-h
tr_da_ts <- ts(da_ts[1:eotr], frequency = freq, start = start(da_ts))
xmin <- time(da_ts)[eotr] - (npts+1) / freq
xmax <- time(da_ts)[eotr] + (h+1) / freq
xlim = c(xmin, xmax)
ylim = c(0, 30)
cat(xmin, ";", xmax)
# Label 1: Actual Earning Line
plot(da_ts, 
     xlim = xlim,
     ylim = ylim,
     main = "JnJ QtlyEarning Multiplicative\nSeasonal Model Forecasting (Fig. 2.15)", 
     xlab = "Year", ylab = "QtlyEarning", lty=1
)
# Plot forecast mean (prepend the last observed data in the training dataset)
dummy_1st_fmean_ts <- ts(
    c(c(da_ts[eotr]), as.numeric(e_qer_arr)), 
    frequency = freq, start = end(tr_da_ts)
)
# Label 2: NULL
lines(dummy_1st_fmean_ts)
# Label 3: Forecast Mean
points(dummy_1st_fmean_ts, pch = 1)
# Plot confidence interval (95%)
t_stats = qnorm(p = 0.975, mean = 0, sd = 1)
dummy_1st_flower_ts <- ts(
    exp(c(c(lg_da_ts[eotr]), e_lg_qer_arr)-t_stats*c(c(0), sqrt(v_qer_arr))), # lower CI for log-normal
    frequency = freq, start = end(tr_da_ts)
)
dummy_1st_fupper_ts <- ts(
    exp(c(c(lg_da_ts[eotr]), e_lg_qer_arr)+t_stats*c(c(0), sqrt(v_qer_arr))), # upper CI for log-normal
    frequency = freq, start = end(tr_da_ts)
)
# Label 4: Forecast 95% Lower Bound 
lines(dummy_1st_flower_ts, lty = 2)
# Label 5: Forecast 95% Upper Bound
lines(dummy_1st_fupper_ts, lty = 2)
# Plot original data
orig_plot_ts <- ts(
    da_ts[(eotr - npts + 1):length(da_ts)], 
    frequency = freq, start = time(da_ts)[eotr] - (npts - 1) / freq
)
# Label 6: NULL
lines(orig_plot_ts)
# Label 7: Actual Earning
points(orig_plot_ts, pch = 19)
legend(
    "topleft", 
    legend = c(
        "Actual Earning Line", NULL, "Forecast Mean", "Forecast 95% Lower Bound",
        "Forecast 95% Upper Bound", NULL, "Actual Earning"
    ), 
#     col = c("black", "red", "blue"), 
    lty = c(1, NA, 2, 2, NA),
    pch = c(NA, 1, NA, NA, 19)
)
1977 ; 1981
No description has been provided for this image

Forecast: Directly apply multiplicative seasonal model

  • The forecast mean from my calculation and the R package are close to each other.
  • This is gives much closer CI to the textbook results than my calculation. My calculation inflate the forecasting variance.
In [144]:
length(jnj_lg_er_ts)
84
In [70]:
npts = 6
h = 8
eotr = length(jnj_lg_er_ts)-h
freq = 4
order = c(0,1,1) # From the solution of 2-6
seasonal = list(order = order, period = freq)
fixed = NULL
include.mean = F
multi_seas_mod_jnj_ts = ts(jnj_lg_er_ts, frequency = freq, start = c(1960, 1))
multi_seas_mod_jnj_fc_res = plot_arima_forecast_fig(
    da_ts=multi_seas_mod_jnj_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', include.mean=include.mean, transform.pars=TRUE,
    main="Forecasts from ARIMA(0,1,1)(0,1,1)[4] with non-zero mean for\n(1-B)(1-B^4)ln(QtlyEarning) = (1-theta*B)(1-Theta*B^4)a_t of JnJ", 
    xlab="Year", ylab="QtlyEarning", ylim=c(2, 3.2)
)
summary(multi_seas_mod_jnj_fc_res)
[1] 71
1977.25 ; 1981
Forecast method: ARIMA(0,1,1)(0,1,1)[4]

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.6559  -0.3492
s.e.   0.1094   0.1104

sigma^2 estimated as 0.008409:  log likelihood = 68.28,  aic = -132.57

Error measures:
                     ME       RMSE        MAE MPE MAPE      MASE       ACF1
Training set 0.01012177 0.08863074 0.07193184 NaN  Inf 0.4398973 0.04998853

Forecasts:
        Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
1979 Q1       2.585637 2.468121 2.703154 2.405912 2.765363
1979 Q2       2.627074 2.502794 2.751353 2.437005 2.817142
1979 Q3       2.611076 2.480383 2.741768 2.411198 2.810953
1979 Q4       2.364514 2.227709 2.501320 2.155288 2.573741
1980 Q1       2.729193 2.549233 2.909154 2.453968 3.004419
1980 Q2       2.770630 2.578687 2.962573 2.477078 3.064181
1980 Q3       2.754632 2.551411 2.957852 2.443833 3.065430
1980 Q4       2.508070 2.294167 2.721974 2.180933 2.835208
No description has been provided for this image
In [61]:
multi_seas_mod_jnj_fc_tb = comb_forecast_res(
    multi_seas_mod_jnj_fc_res, 
    da_ts = multi_seas_mod_jnj_ts,
    eotr=eotr,
    freq=freq
)
multi_seas_mod_jnj_fc_df = as.data.frame(multi_seas_mod_jnj_fc_tb)
multi_seas_mod_jnj_fc_tb
Forecast method: ARIMA(0,1,1)(0,1,1)[4]

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.6559  -0.3492
s.e.   0.1094   0.1104

sigma^2 estimated as 0.008409:  log likelihood = 68.28,  aic = -132.57

Error measures:
                     ME       RMSE        MAE MPE MAPE      MASE       ACF1
Training set 0.01012177 0.08863074 0.07193184 NaN  Inf 0.4398973 0.04998853

Forecasts:
        Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
1979 Q1       2.585637 2.468121 2.703154 2.405912 2.765363
1979 Q2       2.627074 2.502794 2.751353 2.437005 2.817142
1979 Q3       2.611076 2.480383 2.741768 2.411198 2.810953
1979 Q4       2.364514 2.227709 2.501320 2.155288 2.573741
1980 Q1       2.729193 2.549233 2.909154 2.453968 3.004419
1980 Q2       2.770630 2.578687 2.962573 2.477078 3.064181
1980 Q3       2.754632 2.551411 2.957852 2.443833 3.065430
1980 Q4       2.508070 2.294167 2.721974 2.180933 2.835208
A Time Series: 2 × 4
Qtr1Qtr2Qtr3Qtr4
19792.5856372.6270742.6110762.364514
19802.7291932.7706302.7546322.508070
A Time Series: 2 × 4
Qtr1Qtr2Qtr3Qtr4
19790.091698380.096975570.101980050.10675017
19800.140423890.149773900.158573560.16690994
A Time Series: 2 × 4
Qtr1Qtr2Qtr3Qtr4
19792.6419102.5618682.6980002.301585
19802.7850112.6858052.7738382.451867
A Time Series: 8 × 3
ForecastStd. ErrorActual
1979 Q12.5856370.091698382.641910
1979 Q22.6270740.096975572.561868
1979 Q32.6110760.101980052.698000
1979 Q42.3645140.106750172.301585
1980 Q12.7291930.140423892.785011
1980 Q22.7706300.149773902.685805
1980 Q32.7546320.158573562.773838
1980 Q42.5080700.166909942.451867
In [62]:
multi_seas_mod_jnj_fc_df$Forecast
  1. 2.58563736503486
  2. 2.62707369268672
  3. 2.61107553816595
  4. 2.36451442909211
  5. 2.72919337999976
  6. 2.77062970765161
  7. 2.75463155313085
  8. 2.508070444057

As comparison, the following are expectation from my calculation (my calculation is close to the result from package)

$`1`
2.58563748050113
$`2`
2.62707380269453
$`3`
2.61107565060245
$`4`
2.36451453244944
$`5`
2.72919356620702
$`6`
2.77062988840042
$`7`
2.75463173630834
$`8`
2.50807061815533
In [63]:
multi_seas_mod_jnj_fc_df$`Std. Error`^2
  1. 0.00840859360501197
  2. 0.00940426205922513
  3. 0.0103999305134383
  4. 0.0113955989676515
  5. 0.0197188696416149
  6. 0.0224322219144007
  7. 0.0251455741871865
  8. 0.0278589264599723

As comparison, the following are variance from my calculation (my calculation somhow inflate the variance)

$`1`
0.00840857826352835
$`2`
0.0149193615519221
$`3`
0.0214301448403158
$`4`
0.0279409281287095
$`5`
0.0428324475174856
$`6`
0.0624633460377802
$`7`
0.0820942445580749
$`8`
0.101725143078369
In [64]:
par(bg = "white")
da_ts = jnj_er_ts
lg_da_ts = jnj_lg_er_ts
eotr = length(da_ts)-h
tr_da_ts <- ts(da_ts[1:eotr], frequency = freq, start = start(da_ts))
xmin <- time(da_ts)[eotr] - (npts+1) / freq
xmax <- time(da_ts)[eotr] + (h+1) / freq
xlim = c(xmin, xmax)
ylim = c(0, 30)
cat(xmin, ";", xmax)
# Label 1: Actual Earning Line
plot(da_ts, 
     xlim = xlim,
     ylim = ylim,
     main = "JnJ QtlyEarning Multiplicative\nSeasonal Model Forecasting (Fig. 2.15)", 
     xlab = "Year", ylab = "QtlyEarning", lty=1
)
# Plot forecast mean (prepend the last observed data in the training dataset)
multi_seas_e_qer_arr = exp(multi_seas_mod_jnj_fc_df$Forecast+((multi_seas_mod_jnj_fc_df$`Std. Error`)^2)/2)
dummy_1st_fmean_ts <- ts(
    c(c(da_ts[eotr]), as.numeric(multi_seas_e_qer_arr)), 
    frequency = freq, start = end(tr_da_ts)
)
# Label 2: NULL
lines(dummy_1st_fmean_ts)
# Label 3: Forecast Mean
points(dummy_1st_fmean_ts, pch = 1)
# Plot confidence interval (95%)
t_stats = qnorm(p = 0.975, mean = 0, sd = 1)
multi_seas_mod_jnj_fc_res$lower[,2]
dummy_1st_flower_ts <- ts(
    exp(c(c(lg_da_ts[eotr]), multi_seas_mod_jnj_fc_res$lower[,2])), # lower CI for log-normal
    frequency = freq, start = end(tr_da_ts)
)
dummy_1st_fupper_ts <- ts(
    exp(c(c(lg_da_ts[eotr]), multi_seas_mod_jnj_fc_res$upper[,2])), # upper CI for log-normal
    frequency = freq, start = end(tr_da_ts)
)
# Label 4: Forecast 95% Lower Bound 
lines(dummy_1st_flower_ts, lty = 2)
# Label 5: Forecast 95% Upper Bound
lines(dummy_1st_fupper_ts, lty = 2)
# Plot original data
orig_plot_ts <- ts(
    da_ts[(eotr - npts + 1):length(da_ts)], 
    frequency = freq, start = time(da_ts)[eotr] - (npts - 1) / freq
)
# Label 6: NULL
lines(orig_plot_ts)
# Label 7: Actual Earning
points(orig_plot_ts, pch = 19)
legend(
    "topleft", 
    legend = c(
        "Actual Earning Line", NULL, "Forecast Mean", "Forecast 95% Lower Bound",
        "Forecast 95% Upper Bound", NULL, "Actual Earning"
    ), 
#     col = c("black", "red", "blue"), 
    lty = c(1, NA, 2, 2, NA),
    pch = c(NA, 1, NA, NA, 19)
)
1977 ; 1981
A Time Series: 2 × 4
Qtr1Qtr2Qtr3Qtr4
19792.4059122.4370052.4111982.155288
19802.4539682.4770782.4438332.180933
No description has been provided for this image
In [ ]:

CRSP Decile 1 Index

In [389]:
da = read.table("../AFTS_data/Ch02/m-deciles08.txt", header = T)
da[1:5,]
A data.frame: 5 × 5
dateCAP1RETCAP2RETCAP9RETCAP10RET
<int><dbl><dbl><dbl><dbl>
119700130 0.054383-0.004338-0.073082-0.076874
219700227 0.020264 0.020155 0.064185 0.059512
319700331-0.031790-0.028090-0.004034-0.001327
419700430-0.184775-0.193004-0.115825-0.091112
519700529-0.088189-0.085342-0.085565-0.053193
In [401]:
crsp_decile_1_rtn = da$CAP1RET
crsp_decile_1_rtn_ts = ts(crsp_decile_1_rtn, frequency = 12, start = c(1970, 1))
In [393]:
jan = rep(c(1,rep(0,11)),39) # Create Jan dummy
m1 = lm(crsp_decile_1_rtn~jan)
summary(m1)
Call:
lm(formula = crsp_decile_1_rtn ~ jan)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.30861 -0.03475 -0.00176  0.03254  0.40671 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.002864   0.003333   0.859    0.391    
jan         0.125251   0.011546  10.848   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06904 on 466 degrees of freedom
Multiple R-squared:  0.2016,	Adjusted R-squared:  0.1999 
F-statistic: 117.7 on 1 and 466 DF,  p-value: < 2.2e-16
In [408]:
plot_time_fig(crsp_decile_1_rtn_ts, main = "Fig. 2.16(a)", xlab = "Year", ylab = "Simple Return")
plot_acf(da = crsp_decile_1_rtn, lag.max = 40, main = "Fig. 2.16(b)")
crsp_decile_1_adj_rtn = ts(m1$residuals, frequency = 12, start = c(1970, 1))
plot_time_fig(crsp_decile_1_adj_rtn, main = "Fig. 2.16(c)", xlab = "Year", ylab = "Adjusted Return")
plot_acf(da = m1$residuals, lag.max = 40, main = "Fig. 2.16(d)")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [394]:
m2 = arima(crsp_decile_1_rtn, order=c(1,0,0), seasonal=list(order=c(1,0,1), period=12))
summary(m2)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = crsp_decile_1_rtn, order = c(1, 0, 0), seasonal = list(order = c(1, 
    0, 1), period = 12))

Coefficients:
         ar1    sar1     sma1  intercept
      0.1769  0.9882  -0.9144     0.0118
s.e.  0.0456  0.0093   0.0335     0.0129

sigma^2 estimated as 0.004717:  log likelihood = 584.07,  aic = -1160.14

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
In [396]:
par(bg = 'white')
tsdiag(m2,gof=36)
No description has been provided for this image
In [397]:
m3 = arima(crsp_decile_1_rtn, order=c(1,0,0), seasonal=list(order=c(1,0,1), period=12), include.mean = F)
summary(m3)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = crsp_decile_1_rtn, order = c(1, 0, 0), seasonal = list(order = c(1, 
    0, 1), period = 12), include.mean = F)

Coefficients:
         ar1    sar1     sma1
      0.1787  0.9886  -0.9127
s.e.  0.0456  0.0089   0.0335

sigma^2 estimated as 0.00472:  log likelihood = 583.68,  aic = -1161.36

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

Sec. 2.9

Weekly 1yr and 3yr interest rate

In [20]:
da1 = read.table("../AFTS_data/Ch02/w-gs1yr.txt", header = T)
da3 = read.table("../AFTS_data/Ch02/w-gs3yr.txt", header = T)
da1[1:5,]; da3[1:5,]
A data.frame: 5 × 4
yearmondayrate
<int><int><int><dbl>
119621 53.24
219621123.32
319621193.29
419621263.26
519622 23.29
A data.frame: 5 × 4
yearmondayrate
<int><int><int><dbl>
119621 53.70
219621123.75
319621193.80
419621263.77
519622 23.80

Exploration

In [21]:
r1 = da1$rate
r3 = da3$rate
freq = 52
r1_ts = ts(r1, frequency = freq, start = c(1962, 1))
r3_ts = ts(r3, frequency = freq, start = c(1962, 1))
c1 = diff(r1)
c3 = diff(r3)
c1_ts = ts(c1, frequency = freq, start = c(1962, 2))
c3_ts = ts(c3, frequency = freq, start = c(1962, 2))
In [22]:
par(bg = 'white')
plot(r1_ts, main = "Weekly 1yr or 3yr U.S. Interest Rate (Fig. 2.17)", xlab = "Year", ylab = "Percent")
lines(r3_ts, lty = 2)
legend(
    "topleft",
    legend = c("1yr rate", "3yr rate"),
    lty = c(1, 2)
)

plot(c1_ts, main = "Weekly 1yr U.S. Interest Rate Changes (Fig. 2.20(a))", xlab = "Year", ylab = "Percent")
plot(c3_ts, main = "Weekly 3yr U.S. Interest Rate Changes (Fig. 2.20(b))", xlab = "Year", ylab = "Percent")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [23]:
par(bg = "white")
plot(x = r1, y = r3, main = "Fig. 2.18(a)", xlab = "1yr rate", ylab = "3yr rate")
plot(x = c1, y = c3, main = "Fig. 2.18(b)", xlab = "diff(1yr rate)", ylab = "diff(3yr rate)")
No description has been provided for this image
No description has been provided for this image

Fit models

  • Strong serial correlation in ACF of residuals, indicating that there is a unit-root process.
  • Also conducted ADF unit-root tests. The results show that null-hypothesis cannot be rejected.
    • Null-hypothesis is that the time series has unit-root.
In [5]:
m1 = lm(r3~r1)
summary(m1)
Call:
lm(formula = r3 ~ r1)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.82319 -0.37691 -0.01462  0.38661  1.35679 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.83214    0.02417   34.43   <2e-16 ***
r1           0.92955    0.00357  260.40   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5228 on 2465 degrees of freedom
Multiple R-squared:  0.9649,	Adjusted R-squared:  0.9649 
F-statistic: 6.781e+04 on 1 and 2465 DF,  p-value: < 2.2e-16
In [46]:
par(bg = 'white')
plot(m1$residuals, type = 'l', main = "Simple Linear Model (Fig. 2.19(a))", xlab = "Year")
plot_acf(da = m1$residuals, lag.max = 36, main = "Residuals ACF Plot (Fig. 2.19(b))")
pacf(x = m1$residuals, lag.max = 36, main = "Residuals PACF Plot")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [47]:
help(adfTest)
In [53]:
# "lags=12" is chosen from PACF plot
adf_test_res_1 = adfTest(r1, lags = 12, type = c("ct"))
adf_test_res_3 = adfTest(r3, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_1; adf_test_res_3
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -2.4804
  P VALUE:
    0.3749 

Description:
 Tue Feb 27 04:13:44 2024 by user: 
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -2.0718
  P VALUE:
    0.5479 

Description:
 Tue Feb 27 04:13:44 2024 by user: 
In [30]:
box_test_1 = Box.test(m1$residuals, lag = 12, type = 'Ljung')
box_test_1
	Box-Ljung test

data:  m1$residuals
X-squared = 24787, df = 12, p-value < 2.2e-16
  • The serial correlation is largely reduced, but there are still some serial correlation in ACF of residuals, indicating that the model is not adequate. We need to model some serial correlation in residuals. The ACF plot shows strong serial correlation at lag-1, so that we use MA(1) to model the error.
In [24]:
m2 = lm(c3~-1+c1) # No intercept
summary(m2)
Call:
lm(formula = c3 ~ -1 + c1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42469 -0.03589 -0.00127  0.03456  0.48911 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
c1 0.791935   0.007337   107.9   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06896 on 2465 degrees of freedom
Multiple R-squared:  0.8253,	Adjusted R-squared:  0.8253 
F-statistic: 1.165e+04 on 1 and 2465 DF,  p-value: < 2.2e-16
In [25]:
plot_acf(da = m2$residuals, lag.max = 36, main = "Fig. 2.21(b)")
No description has been provided for this image
In [28]:
# "lags=12" is chosen from PACF plot
adf_test_res_c1 = adfTest(c1, lags = 12, type = c("ct"))
adf_test_res_c3 = adfTest(c3, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_c1; adf_test_res_c3
Warning message in adfTest(c1, lags = 12, type = c("ct")):
“p-value smaller than printed p-value”
Warning message in adfTest(c3, lags = 12, type = c("ct")):
“p-value smaller than printed p-value”
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -13.6578
  P VALUE:
    0.01 

Description:
 Thu Feb 29 12:11:17 2024 by user: 
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -13.933
  P VALUE:
    0.01 

Description:
 Thu Feb 29 12:11:17 2024 by user: 
In [29]:
box_test_2 = Box.test(m2$residuals, lag = 12, type = 'Ljung')
box_test_2
	Box-Ljung test

data:  m2$residuals
X-squared = 141.84, df = 12, p-value < 2.2e-16
  • The serial correlations in ACF plots no longer exists. But the Box-Ljung test is still significant. Maybe we can improve the model more to reduce some serial correlation in the residuals.
In [31]:
m3 = arima(c3, order = c(0,0,1), xreg = c1, include.mean = F)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(m3$residuals^2)/sum(c3^2)
summary(m3); rsq
Call:
arima(x = c3, order = c(0, 0, 1), xreg = c1, include.mean = F)

Coefficients:
         ma1      c1
      0.1823  0.7936
s.e.  0.0196  0.0075

sigma^2 estimated as 0.0046:  log likelihood = 3136.62,  aic = -6267.23

Training set error measures:
                        ME       RMSE        MAE MPE MAPE      MASE
Training set -7.953792e-05 0.06782056 0.04843167 NaN  Inf 0.3766059
                     ACF1
Training set 0.0008170501
0.831007651496559
In [32]:
plot_acf(da = m3$residuals, lag.max = 36, main = "ACF of residuals after taking difference\nand model error as a MA(1) process")
No description has been provided for this image
In [33]:
log(length(c1))
7.81035268372429
In [35]:
box_test_3 = Box.test(m3$residuals, lag = 12, type = 'Ljung')
box_test_3
	Box-Ljung test

data:  m3$residuals
X-squared = 54.422, df = 12, p-value = 2.297e-07

Sec. 2.10

Example of consistent estimation for covariance matrix (1yr and 3yr rates)

  • For heteroskedasticity-consistent or heteroskedasticity&autocorrelation-consistent estimators.
In [32]:
da1 = read.table("../AFTS_data/Ch02/w-gs1yr.txt", header = T)
da3 = read.table("../AFTS_data/Ch02/w-gs3yr.txt", header = T)
da1[1:5,]; da3[1:5,]
A data.frame: 5 × 4
yearmondayrate
<int><int><int><dbl>
119621 53.24
219621123.32
319621193.29
419621263.26
519622 23.29
A data.frame: 5 × 4
yearmondayrate
<int><int><int><dbl>
119621 53.70
219621123.75
319621193.80
419621263.77
519622 23.80
In [33]:
r1 = da1$rate
r3 = da3$rate
freq = 52
r1_ts = ts(r1, frequency = freq, start = c(1962, 1))
r3_ts = ts(r3, frequency = freq, start = c(1962, 1))
c1 = diff(r1)
c3 = diff(r3)
c1_ts = ts(c1, frequency = freq, start = c(1962, 2))
c3_ts = ts(c3, frequency = freq, start = c(1962, 2))
In [34]:
m1 <- lm(c3~c1)
summary(m1)
Call:
lm(formula = c3 ~ c1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42459 -0.03578 -0.00117  0.03467  0.48921 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.0001051  0.0013890  -0.076     0.94    
c1           0.7919323  0.0073391 107.906   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06897 on 2464 degrees of freedom
Multiple R-squared:  0.8253,	Adjusted R-squared:  0.8253 
F-statistic: 1.164e+04 on 1 and 2464 DF,  p-value: < 2.2e-16
  • HC estimator
    • GPT: R use white correction in OLS summary
In [37]:
help(vcovHC)
In [38]:
# Applying White's correction
# `type` parameter 'arg' should be one of “HC3”, “const”, “HC”, “HC0”, “HC1”, “HC2”, “HC4”, “HC4m”, “HC5”
hc1_robust_se <- vcovHC(m1, type = "HC1") # HC1 is often recommended; it's a variant of White's correction
hc1_robust_se
A matrix: 2 × 2 of type dbl
(Intercept)c1
(Intercept) 1.928016e-06-4.212675e-07
c1-4.212675e-07 2.672754e-04
In [43]:
help(coeftest)
In [42]:
# Getting the summary with corrected standard errors
summary_with_hc1_robust_se <- coeftest(m1, hc1_robust_se)
coeftest(m1, vcovHC);
coeftest(m1, vcovHC, type = "HC0");
print(summary_with_hc1_robust_se);
dwtest(m1);
jarque.bera.test(m1$residuals);
Box.test(m1$residuals, lag = 33, type = 'Ljung')
t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00139139 -0.0756   0.9398    
c1           0.79193231  0.01659956 47.7080   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00138797 -0.0758   0.9396    
c1           0.79193231  0.01634193 48.4601   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00138853 -0.0757   0.9396    
c1           0.79193231  0.01634856 48.4405   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

	Durbin-Watson test

data:  m1
DW = 1.6456, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0
	Jarque Bera Test

data:  m1$residuals
X-squared = 1644.6, df = 2, p-value < 2.2e-16
	Box-Ljung test

data:  m1$residuals
X-squared = 230.05, df = 33, p-value < 2.2e-16
  • HAC estimator
    • GPT: Use heteroscedasticity and autocorrelation consistent (HAC) estimator in R
In [45]:
help(vcovHAC)
In [46]:
# Compute HAC standard errors
hac_se <- vcovHAC(m1)
hac_se
A matrix: 2 × 2 of type dbl
(Intercept)c1
(Intercept) 2.505407e-06-1.306652e-06
c1-1.306652e-06 3.051041e-04
In [47]:
# Summary of the model with HAC standard errors
# Doesn't have effect
summary(m1, robust = hac_se)
Call:
lm(formula = c3 ~ c1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42459 -0.03578 -0.00117  0.03467  0.48921 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.0001051  0.0013890  -0.076     0.94    
c1           0.7919323  0.0073391 107.906   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06897 on 2464 degrees of freedom
Multiple R-squared:  0.8253,	Adjusted R-squared:  0.8253 
F-statistic: 1.164e+04 on 1 and 2464 DF,  p-value: < 2.2e-16
In [48]:
# Alternatively, you can directly use coeftest() from lmtest to test coefficients with HAC standard errors
coeftest(m1, vcov. = hac_se)
t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00158285 -0.0664    0.947    
c1           0.79193231  0.01746723 45.3382   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  • The following function gives closer answer at the pp100 in the textbook.
In [49]:
coeftest(m1, NeweyWest(m1, lag = 12, prewhite = F)) # From the solution to 2-7
t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00149696 -0.0702    0.944    
c1           0.79193231  0.02031114 38.9901   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  • Including lag in linear regression
    • GPT: R linear model using lag operator
In [83]:
help(merge)
In [93]:
# Create a zoo object for the time series
c3_zoo <- zoo(c3)
c1_zoo <- zoo(c1)

# Lag the series by 1 period
c3_lag1_zoo <- lag(c3_zoo, -1, na.pad=TRUE)
c1_lag1_zoo <- lag(c1_zoo, -1, na.pad=TRUE)

# Combine the original and lagged series
c3_c1_zoo <- merge(c3_zoo, c1_zoo, c3_lag1_zoo, c1_lag1_zoo, all = TRUE)
dim(c3_c1_zoo); c3_c1_zoo[1:5,]
  1. 2466
  2. 4
  c3_zoo c1_zoo c3_lag1_zoo c1_lag1_zoo
1   0.05   0.08          NA          NA
2   0.05  -0.03        0.05        0.08
3  -0.03  -0.03        0.05       -0.03
4   0.03   0.03       -0.03       -0.03
5  -0.04   0.00        0.03        0.03
In [111]:
help(lm)
In [113]:
# Fit a linear model
model_zoo <- lm(c3_zoo ~ c1_zoo+c3_lag1_zoo+c1_lag1_zoo, data = c3_c1_zoo, na.action = na.exclude)
# View the summary of the model
summary(model_zoo)
dwtest(model_zoo);
jarque.bera.test(model_zoo$residuals);
Box.test(model_zoo$residuals, lag = 33, type = 'Ljung')
Call:
lm(formula = c3_zoo ~ c1_zoo + c3_lag1_zoo + c1_lag1_zoo, data = c3_c1_zoo, 
    na.action = na.exclude)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.44812 -0.03551 -0.00075  0.03408  0.45821 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.686e-05  1.367e-03  -0.064    0.949    
c1_zoo       7.971e-01  7.692e-03 103.632   <2e-16 ***
c3_lag1_zoo  1.766e-01  1.983e-02   8.906   <2e-16 ***
c1_lag1_zoo -1.580e-01  1.745e-02  -9.058   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06785 on 2461 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.8312,	Adjusted R-squared:  0.831 
F-statistic:  4039 on 3 and 2461 DF,  p-value: < 2.2e-16
	Durbin-Watson test

data:  model_zoo
DW = 1.9865, p-value = 0.3649
alternative hypothesis: true autocorrelation is greater than 0
	Jarque Bera Test

data:  model_zoo$residuals
X-squared = 1620.5, df = 2, p-value < 2.2e-16
	Box-Ljung test

data:  model_zoo$residuals
X-squared = 131.6, df = 33, p-value = 1.009e-13

Prove that the equation on pp101 is the $j^{th}$ diagnoal element of the equation (2.49).

  • From equation (2.49) $$ \begin{eqnarray} Cov(\hat{\boldsymbol{\beta}})_{HC} = \left[\sum_{t=1}^T \mathbf{x}_t\mathbf{x}_t'\right]^{-1} \left[\sum_{t=1}^T \hat{e}_t^2\mathbf{x}_t\mathbf{x}_t'\right] \left[\sum_{t=1}^T \mathbf{x}_t\mathbf{x}_t'\right]^{-1} \end{eqnarray} $$

  • The equation on pp101 $$ Var(\hat{\beta}_j)_{HC} = \frac{\sum_{t=1}^T\hat{e}_t^2\hat{v}_t^2}{\left(\sum_{t=1}^T\hat{v}_t^2\right)^2} $$

where $\hat{v}_t$ is the least-squares residual the auxiliary regression $x_{jt}=\mathbf{x}_{-j,t}'\boldsymbol{\gamma}+v_t, t=1,\cdots,T$.

  • Basically, we want to prove that $$ \begin{eqnarray} \mathbf{e}_j'\left[\sum_{t=1}^T \mathbf{x}_t\mathbf{x}_t'\right]^{-1} \left[\sum_{t=1}^T \hat{e}_t^2\mathbf{x}_t\mathbf{x}_t'\right] \left[\sum_{t=1}^T \mathbf{x}_t\mathbf{x}_t'\right]^{-1}\mathbf{e}_j = \frac{\sum_{t=1}^T\hat{e}_t^2\hat{v}_t^2}{\left(\sum_{t=1}^T\hat{v}_t^2\right)^2} \end{eqnarray} $$

where $\mathbf{e}_j$ is the unit vector with $j^{th}$ element being $1$ and other element being $0$.

Before proving the general case, we can first check if this is true when all the columns of $\mathbf{X}$ are orthogonal.

We can write $\left[\sum_{t=1}^T\mathbf{x}_t\mathbf{x}_t'\right]=\mathbf{X}'\mathbf{X}$ as a block matrix and use the formula of inversion of a block matrix to prove this.

Sec. 2.11

Daily simple return of value-&equal-weighted index from CRSP

In [117]:
da_daily = read.table("../AFTS_data/Ch02/d-ibm3dx7008.txt", header = T)
vwretd = da_daily$vwretd
ewretd = da_daily$ewretd
In [135]:
par(bg = "white")
plot(acf(abs(vwretd), lag.max = 300, plot = F), main = "ACF of abs(vwretd) (Fig. 2.22(a))")
plot(acf(abs(ewretd), lag.max = 300, plot = F), main = "ACF of abs(ewretd) (Fig. 2.22(b))")
No description has been provided for this image
No description has been provided for this image
In [ ]: